Thermal and nuclear quantum effects in the hydrogen bond dynamical symmetrization phase transition of AlOOH
Abstract
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 OH 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 OO 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 HO 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 highsymmetry and protoncentered phase, which is brought about by the onset of proton tunneling, could be confirmed by quasielastic 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 highpressure 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 (SanoFurukawa 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 Xray 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 posttransition 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 secondorder phase transition from an asymmetric (hydrogen offcentered  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 VIIX 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; SanoFurukawa et al., 2008; Tsuchiya et al., 2008); moreover, AlOOH belongs to the spacegroup, 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 doublewell picture stands for AlOOH as well. In particular, contrarily to the highpressure ice structure, the oxygenoxygen midpoint is not a symmetry center in AlOOH, whereas it is in the spacegroup; 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 Hbonds 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 Xray 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 Langevintype 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:
(1) 
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):
(2)  
(3) 
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 (DFTbased) 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; BedoyaMartinez et al., 2014; HernandezRojas 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 pathintegralbased 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 OH 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 BeckeLeeYangParr functional (Becke, 1988; Lee et al., 1988). The interaction between the ionic cores and the valence electrons is described through ultrasoft pseudopotentials with nonlinear core corrections. The KohnSham orbitals are expanded in planewaves with Ry energy cutoff. Our simulation cell is a supercell containing AlOOH units. The corresponding Brillouin Zone is sampled by a MonkhorstPack 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 Xray 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 
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 Xray 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): offcentered 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 abinitio 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 abinitio 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 Xray experiments showed a change of compressibility of AlOOH around GPa (SanoFurukawa 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 lowpressure (LP) structure, the other for the highpressure (HP) structure (see Figure 2), as done in previous studies (Tsuchiya et al., 2002; SanoFurukawa 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 (SanoFurukawa 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 OHO 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 OO average distance is a relevant parameter for any transition from asymmetric OHO to symmetric OHO bonds; indeed, in pure ice, proton tunneling and the onset of symmetrization happen at OO distances around Å (Benoit & Marx, 2005), and this is generally accepted as a good criterion for the onset of symmetrization. The average OO 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 zeropoint motion. Our results are in agreement with experimental data (Komatsu et al., 2006; SanoFurukawa et al., 2008; Kuribayashi et al., 2014), though we slightly overestimate the OO distance as pressure increases to about GPa with respect to very recent Xray data (Kuribayashi et al., 2014). One may note that both experimental results and ours for the OO distance, remain significantly above the 2.4 Å threshold until over 25 GPa.
In order to investigate the hydrogen bonds in AlOOH, we compute the oxygenhydrogen (OH) 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 OH bonds and around Å for hydrogen HO bonds. These distances are in agreement with experimental measurements on AlOOD obtained through neutron scattering around GPa (SanoFurukawa 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 OO mean distance decreases (see Figure 3); this is correlated with a dilatation of the OH 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 oxygenoxygen distance is larger than 2.4Å (Figure 3), while the OH distance at the same pressure is less than 1.1Å, that is, less than half the OO 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 OH distances between a given proton H and its two nearestneighbor 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 OO segment.
Figure 5 shows the distribution of the proton transfer coordinate computed in two different ways: from the timeaveraged atomic positions and from the instantaneous atomic positions. The timeaveraged 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 ionocovalent OH 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 (SanoFurukawa 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 protonordered asymmetric hydrogen bond to a dynamically disordered symmetric hydrogen bond taking place around 1520 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 protondisordered asymmetric hydrogen bond to a protondisordered 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 Xray scattering, in which proton disorder should be barely visible, only the averaged symmetry can be seen at high pressure. Proton tunneling at the timescale of the experiment thus triggers the dynamical symmetrization of the hydrogen bond that leads to the hydrogencentered (HC) phase, which is observed experimentally (SanoFurukawa 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 Xray diffraction experiments, and paircorrelation 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 timedependent 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 highfrequency OH stretching modes (around cm at GPa), and four OHO bending modes between and cm. The two lowfrequency bending modes vibrate in the direction, while the two bending modes at higher frequencies are in the plane, which contains the OHO bonds. When pressure increases, the OH 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 OH bending modes (around cm at ) are in the plane and the two lowest ones along vibrate in the  cm range. In contrast, the OH 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 secondorder transition from the HOC to the HC structure in AlOOH, which is accompanied by a strong softening of the OH 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 OH 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 OO direction, while is orthogonal to in the plane. In the configuration, the OH stretching mode eigenvector forms a small angle with the OO 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 modespecific 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:
(4) 
where is the shortest OH distance (corresponding to the OH 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:
(5) 
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 OH stretching modes and the OHO 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 highfrequency OH 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 OH pair correlation function. Upon compression, the OH stretching modes soften gradually, their frequencies dropping to approximately and cm. Moreover, the width of the highfrequency 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 OH stretching mode frequency ( cm) does not vary appreciably with pressure, at least up to GPa, while the bending mode frequencies increase continuously. The second OH 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 OH stretching mode up to about GPa, where bending and stretching modes mix up, is consistent with a secondorder phase transition as also guessed by the discontinuity of the compressibility. This is also quite close to the highpressure behavior of the OH 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 OHO 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 protondisordered 2nm phase with asymmetric hydrogen bonds to a protondisordered phase with symmetric hydrogen bonds. The proton disorder decreases with pressure and a protonordered 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 highpressure 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 semiclassical atomic trajectories fulfilling the BoseEinstein 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 lowpressure hydrogen offcentered (HOC) phase with space group to a more symmetric, hydrogencentered (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 protonordered HC phase. The apparent contradiction with recent Xray diffraction data (Kuribayashi et al., 2014) providing experimental evidence that the transition to the highsymmetry 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 OO distances down to Å  the usual threshold for the hydrogenbond symmetrization (Benoit & Marx, 2005; Morrone et al., 2009; Lin et al., 2011). Proton tunneling implies that the timeaveraged proton site coincides with the OO 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 DebyeWaller factors through quasielastic neutron scattering (Teixeira J. & Luzar A., 1999) could distinguish between the two regimes.
As in other orderdisorder phase transitions, the analysis of the OH modes is crucial to point out the transition. From the trends of Hrelated 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 OH stretching modes considerably weaken and eventually fade out. Those trends are fully consistent with the analysis of the structural observables (average positions and paircorrelation 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 highresolution 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.
Acknowledgments:
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’ÎledeFrance (DIM Oxymore).
This work was granted access to the HPC resources of CINES under the allocation 2016096719 made by GENCI.
References
 Aoki et al. (1996) Aoki, K., Yamawaki, H., Sakashita, M. & Fujihisa, H. (1996): Infrared absorption study of the hydrogenbond symmetrization in ice to 110 GPa, Phys. Rev. B 54, 1567315677
 Baroni et al. (2001) Baroni, S., de Gironcoli, S., Dal Corso, A. & Giannozzi, P. (2001): Phonons and related crystal properties from densityfunctional perturbation theory, Rev. Modern Phys. 73, 515562
 Basire et al. (2013) 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, 1259112601
 Becke (1988) Becke, A. D. (1988): Densityfunctional exchangeenergy approximation with correct asymptotic behavior, Phys. Rev. A 38, 30983100
 BedoyaMartinez et al. (2014) BedoyaMartinez, 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
 Benoit et al. (1998) Benoit, M., Marx, D. & Parrinello, M. (1998): Tunnelling and zeropoint motion in highpressure ice, Nature 392, 258261
 Benoit & Marx (2005) Benoit M. & Marx, D. (2005): The shapes of protons in hydrogen bonds depend on the bond length, Chem. Phys. Phys. Chem. 6, 17381741
 Berne & Thirumalai (1986) Berne, B. J.& Thirumalai, D. (1986): On the simulation of quantum systems: path integral methods, Ann. Rev. Phys. Chem. 37, 401424
 Brieuc et al. (2016) Brieuc, F., Bronstein, Y., Dammak, H., Depondt, P., Finocchi, F. & Hayoun, M. (2016): Zeropoint Energy Leakage in Quantum Thermal Bath Molecular Dynamics Simulations, J. Chem. Theory and Comput. 12, 56885697
 Bronstein et al. (2014) Bronstein, Y., Depondt, P., Finocchi, F. & Saitta A. M. (2014): Quantumdriven phase transition in ice described via an efficient Langevin approach, Phys. Rev. B 89, 214101
 Bronstein et al. (2016) 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
 Callen & Welton (1951) Callen, H. B. & Welton, T. A. (1951): Irreversibility and generalized noise, Phys. Rev. 83, 3440
 Calvo et al. (2012) Calvo, F., VanOanh, N. T., Parneix, P. & Falvo, C. (2012): Vibrational spectra of polyatomic molecules assisted by quantum thermal baths, Phys. Chem. Chem. Phys. 14, 1050310506
 Dakhnovskii et al. (1995) 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
 Dammak et al. (2009) Dammak, H., Chalopin, Y., Laroche, M., Hayoun, M. & Greffet, J. J. (2009): Quantum thermal bath for molecular dynamics simulation, Phys. Rev. Lett. 103, 190601
 Dammak et al. (2012) 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
 Friedrich et al. (2007) Friedrich, A., Haussühl, E., Boehler, R., Morgenroth, W., JuarezArellano, E. A. & Winkler, B. (2007): Singlecrystal structure refinement of diaspore at 50 GPa, Am. Mineral. 92, 16401644
 Giannozzi et al. (2009) Giannozzi, P. et al. (2009): QUANTUM ESPRESSO: a modular and opensource software project for quantum simulations of materials, J. Phys.: Cond. Matter 21, 395502
 Goncharov et al. (1996) 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 hydrogenbonded phase, Science 273, 218220
 Goncharov et al. (1999) 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 19982001
 Hermann & Mookerjee (2016) Hermann, H. & Mookerjee, M. (2016): Highpressure phase of brucite stable at earth’s mantle transition zone and lower mantle conditions, PNAS DOI: 10.1073/pnas/1611571113
 HernandezRojas et al. (2015) HernandezRojas, J., Calvo, F. & Noya, E. G.: (2015) Applicability of Quantum Thermal Baths to Complex ManyBody Systems with Various Degrees of Anharmonicity, J. Chem. Theory Comput. 11, 861870
 Hills (1979) Hills, R. J. (1979): Crystal structure refinement and electron density distribution in diaspore, Phys. Chem. Minerals 5, 179200
 Izaguirre et al. (2001) Izaguirre, J. A., Catarello, D. P., Wozniak, J. M. & Skeel, R. D. (2001): Langevin stabilization of molecular dynamics, J. Chem. Phys. 114, 2090
 Jahn et al. (2012) Jahn, S., Wunder, B., KochMüller, M., Tarrieu, L., Pöhle, M., Watenphul, A. & Taran, M. N. (2012): Pressureinduced hydrogen bond symmetrisation in guyanaite, Î²CrOOH: evidence from spectroscopy and ab initio simulations, Eur. J. Mineral. 24, 839850
 Johannsen (1998) Johannsen, P. G. (1998): Vibrational states and optical transitions in hydrogen bonds, J. Physics: Cond. Matter 10, 22412260
 Komatsu et al. (2006) Komatsu, K., Kuribayashi, T., Sano, A., Ohtani, E. & Kudoh, Y. (2006): Redetermination of the highpressure modification of AlOOH from singlecrystal synchrotron data, Acta Cryst. E 62, i216i218
 Kubo et al. (1978) Kubo, R., Toda, M. & Hashitsume, N. (1978): Statistical Physics II: Nonequilibrium statistical mechanics, Springer
 Kuribayashi et al. (2014) Kuribayashi, T., SanoFurukawa, A. & Nagase, T. (2014): Observation of pressureinduced phase transition of Î´AlOOH by using singlecrystal synchrotron Xray diffraction method, Phys. Chem. Minerals 41, 303312
 Li et al. (2006) Li, S., Ahuja, R. & Johansson, B. (2006): The elastic and optical properties of the highpressure hydrous phase AlOOH, Solid State Comm. 137, 101105
 Lee et al. (1988) Lee, C., Yang, W. & Parr, R. G. (1988): Development of the ColleSalvetti correlationenergy formula into a functional of the electron density, Phys. Rev. B 37, 785789
 Lin et al. (2011) Lin, L., Morrone, J. A. & Car, R. (2011): Correlated tunneling in hydrogen bonds, J. Stat. Phys. 145, 365384
 Maradudin et al. (1990) 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
 Marshall & Lovesey (1971) Marshall, W. & Lovesey, S. W. (1971): Theory of thermal neutron scattering: the use of neutrons for the investigation of condensed matter. Oxford University Press
 Marx & Parrinello (1996) Marx, D. & Parrinello, M. (1996): Ab initio path integral molecular dynamics: Basic ideas, J. Chem. Phys. 104, 40774082
 Morrone et al. (2009) 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
 Nachimuthu et al. (2012) 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, 812
 Nielsen & Martin (1985) Nielsen, O. H. & Martin, R. M. (1985): Quantummechanical theory of stress and force, Phys. Rev. B 32, 37923805
 Ohira et al. (2014) 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, 1217
 Ohtani et al. (2001) 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, 39913993
 Panero & Stixrude (2004) 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, 421431
 Perdew et al. (1996) Perdew, J. P., Burke, K. & Ernzerhof, M. (1996): Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865
 Sano et al. (2004) Sano, A., E. Ohtani, Kubo, T. & Funakoshi, K. (2004): In situ Xray observation of decomposition of hydrous aluminum silicate AlSiO OH and aluminum oxide hydroxide AlOOH at high pressure and temperature, J. Phys. Chem. Solids 65, 15471554
 Sano et al. (2008) 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 coremantle boundary, Geophys. Research Lett. 35, L03303
 SanoFurukawa et al. (2008) SanoFurukawa, 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, 15581567
 SanoFurukawa et al. (2009) SanoFurukawa, 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 hydrogenbond symmetrization, Am. Mineral. 94, 12551261
 Savin et al. (2012) Savin, A. V., Kosevich, Y. A. & Cantarero, A. (2012): Semiquantum molecular dynamics simulation of thermal properties and heat transport in lowdimensional nanostructures, Phys. Rev. B 86, 064305
 Struzhkin et al. (1997) 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 44464449
 Sugimura et al. (2008) Sugimura, E., Iiataka, T., Hirose, K., Kawamura, K., Sata, N. & Ohishi, Y. (2008): Compression of HO ice and implications for hydrogenbond symmetrization: Synchrotron xray diffraction measurements and densityfunctional calculations, Phys. Rev. B 77 214103
 Suzuki et al. (2000) 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) 689693
 Teixeira J. & Luzar A. (1999) Teixeira J. & Luzar A. (1999): Physics of liquid water in Hydration Processes in Biology: Theoretical and Experimental Approaches, M.C. BellisentFunel ed., IOS Press, NATO Science Studies, 3565
 Tsuchiya et al. (2002) Tsuchiya, J., Tsuchiya, T., Tsuneyuki, T. & Yamanaka, T. (2002): First principles calculation of a highpressure hydrous phase, AlOOH, Geophys. Research Lett. 29, 1909
 Tsuchiya et al. (2008) Tsuchiya, J., Tsuchiya; T. & Wentzcovitch, R. M. (2008): Vibrational properties of AlOOH under pressure, Am. Mineral. 93, 477482
 Tsuchiya & Tsuchiya (2009) Tsuchiya, J. & Tsuchiya, T. (2009): Elastic properties of AlOOH under pressure: First principles investigation, Phys. Earth and Planetary Interiors 174, 122127
 Tsuchiya (2013) 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, 45704573
 Umemoto et al. (2015) 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
 Vanpeteghem et al. (2002) 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, 11191121
 Vinet et al. (1987) 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
 Winkler et al. (2001) Winkler, B., Hytha, M., Pickard, P., Milman, V., Warren, M. & Segall, M. (2001): Theoretical investigation of bonding in diaspore, Eur. J. Mineral. 13, 343349