Quantum driven proton diffusion in brucite-like minerals under high pressure

We investigate the elementary steps at the microscopic level for proton diffusion in brucite under high pressure, which results from a complex interplay between two processes: the O-H reorientations motion around the $\mathbf c$ axis and O-H covalent bond dissociations. First-principle path-integral molecular dynamics simulations reveal that the increasing pressure tends to lock the former motion, while, in contrast, it activates the latter which is mainly triggered by nuclear quantum effects. These two competing effects therefore give rise to a pressure sweet spot for proton diffusion within the mineral. In brucite \ce{Mg(OH)2}, proton diffusion reaches a maximum for pressures close to 70GPa, while the structurally similar portlandite \ce{Ca(OH)2} never shows proton diffusion within the pressure range and time scale that we explored. We analyze the different behaviors of brucite and portlandite, which might constitute two prototypes for other minerals with the same structure.


Introduction
Hydroxide minerals play an important role in several problems in geology, surface science or for industrial applications. Among them, brucite Mg(OH) 2 can be formed at the interface between periclase MgO and water at ambient conditions [1][2][3] . The trigonal brucite structure consists of alternating layers along the c axis that terminate with hydroxyls ( figure 1). This structure is common to other hydroxides of divalent metals, such as Ca(OH) 2 , Ni(OH) 2 and Cd(OH) 2 . Portlandite Ca(OH) 2 is the main component of cements and concretes, which motivated a large number of investigations about its elastic properties. Because of their anisotropic structure, brucite isostructural minerals are much more compressible along the c axis than in the other two directions, parallel to the stacks.
Mg(OH) 2 can also act as a water vector in subduction zones, through complex processes that take place within the Earth interior 4,5 . Therefore, the behavior of brucite and brucite-like minerals at very high pressure has been widely investigated. X-ray diffraction of Mg(OH) 2 up to 78 GPa showed that the c/a ratio decreases steadily from ambient pressure up to about 25 GPa and then remains almost constant 6 . Those results suggest that the properties of brucite at very high pressure, and in particular the nature of the inter-layer bonding, could differ significantly from ambient conditions. Moreover, the hydroxyl groups that are mostly parallel to the c axis at ambient conditions slant in three equivalent positions as the inter-layer distance shrinks even under moderate pressure 7,8 or when decreasing temperature 9 . Besides reducing the global symmetry from P3m1 down to P3, the slanted OH groups can significantly alter several physical properties of brucites. Firstly, it could allow for the formation of hydrogen bonds between the layers 10 , which eventually reinforce under further compression and modify the compressibility along c 11 . Secondly, when the protons form a non null θ angle with the c axis, they cannot arrange in a static ordered structure. Such proton disorder, which is closely related to proton frustration 12,13 , has also been invoked as the reason for the pressure-induced hydrogen sublattice amorphization in brucites 7,14 . The existence of a quasi two-dimensional proton liquid in those extreme conditions can be conjectured, but the properties of the whole structure, if stable, have so far escaped a precise characterization. In particular, the occurrence of proton hopping is plausible, but whether this process results in a long-range diffusion is a totally open question.
From the theoretical viewpoint, the previous observations call for a dynamical treatment of the proton arrangement within the brucite structure at high pressure. Moreover, in such conditions nuclear quantum effects (NQE), that is, all the properties that go beyond a purely classical description of ion dynamics 15 , such as zero-point energy, ion delocalization and tunneling, can have an important impact on the brucite properties. Its light mass makes hydrogen prone to fast diffusion and subject to significant quantum effects even at room temperature.
In particular, the high O-H stretching frequency implies a non-negligible zero-point energy of about 0.2 eV that is crucial when hopping through local sites which are separated by barriers that classical nuclei cannot overcome simply by thermal fluctuations [16][17][18] . Indeed, the seminal paper by Benoit et al. 16 Figure 1. Left-hand side: description of brucite minerals structure. The primitive cell is shaded; the simulation box is defined in-plane by vectors a = a(3/2, √ 3/2, 0) and b = a(0, √ 3, 0), and contains 3 Mg(OH) 2 units. The hydrogen atoms in the P3m1 symmetry structure (light grey) are in the 2d Wyckoff sites, while the hydrogen atoms in the P3 structure (dark grey) are in the disordered 6i Wyckoff sites. The latter is the stable structure at the pressures of our simulations. Right-hand side: Proton diffusion mechanism which requires a succession of two alternate steps: proton reorientation (green arrows) i.e. O-H rotation around c; and dissociation (blue arrows) i.e. proton hopping between oxygen layers. Oxygen atoms are colored in red, hydrogen atoms in grey. The size of atoms decreases with depth.  Table 1. Lattice constants of brucite from our simulations and c/a ratios. The numerical errors are of approximately 0.01Å for a, 0.02Å for c and therefore of the order of 10 −2 for c/a. becomes comparable to that along a 45 , in connection with the formation of inter-layer hydrogen bonds. This effect is expected to enhance the cohesion of the system and thus to stabilize the structure, which might be connected to the existence of a rather large barrier hindering brucite decomposition. Moreover, the formation of inter-layer hydrogen bonds could allow the protons to hop between the two oxygen atoms on the facing Mg(OH) 2 layers. All along the paper, we analyze in detail the reason behind the increasing inter-layer repulsion and discuss the role of the nuclear quantum effects in the formation of hydrogen bonds between the layers.

Proton diffusion mechanism
As in other hydrous minerals 19,20 , hydrates 46 and ices 16,47 the proton can hop along hydrogen bonds. Following Dupuis et al. 33 , we refer to this process here as dissociation, as it implies the breaking of a covalent O-H bond to form another distinct O-H covalent bond. However, another process, namely the reorientation (the hopping between the three 6i Wyckoff sites in the P3 space group), is necessary for the proton to move away from the initial O site through the crystal. Only the combination of dissociation and reorientation can drive proton diffusion; one of the two processes alone would simply imply a back-and-forth proton motion between neighboring sites. Proton diffusion in brucite-like minerals is therefore a two-step compound process, intrinsically different from the Grotthuss mechanism that is active in water and other H-bonded hydrates 37,48 .
In P3 Mg(OH) 2 , the proton experiences an effective triple well potential among the 6i Wyckoff sites, which controls reorientation within the (ab) plane. We refer to these reorientation events as "in-plane" motion. On the other hand, the "out-of-plane" dissociation mechanism involves an effective double well potential along the O-O direction characterizing the covalent and hydrogen bonds. A recent study 33 analyzed the proton diffusion mechanism in portlandite for different temperatures and concluded that thermal activation of the reorientation hopping can be the limiting factor for proton diffusion. Although the authors guessed that nuclear quantum effects could favor the dissociation mechanism, the latter process was studied only classically.

Results
In the following, we mainly focus on the effect of pressure in brucite and unravel the complex and quantum driven proton diffusion in this material, accounting for NQE in a consistent framework for both reorientation and dissociation processes. Rising pressure tends to increase the strength of hydrogen bonds. As a consequence, OH dissociation and proton reorientation display opposite trends with pressure. Furthermore, we show that dissociation is significantly enhanced by NQE while reorientation is mostly thermally activated. Finally, the behavior of protons in portlandite is analyzed, and differences between Mg(OH) 2 and Ca(OH) 2 are discussed.

Quasi-2D proton layer
Before addressing proton diffusion, we need to examine the structural changes of the hydrogen layers under increasing pressure and the possible formation of an almost two-dimensional proton layer. Figure 2 shows the probability distribution of protons along the c direction. Initially, each proton lies either in the upper or lower plane -we therefore label the protons according to their initial location and compute the corresponding "lower-layer" or the "upper-layer" proton distributions. As pressure increases, the overall distribution width decreases, the double-peak structure still visible at 30 GPa disappears as the lower-and the upper-layer distributions tend to merge (see Fig. 2). Accounting for NQE via PIMD is crucial for the study of this process: zero-point energy yields an important contribution to the proton distribution width and makes the merging of the two layers occur at lower pressure than for classical protons.
An other interesting effect appears when considering hydrogen motion between the two layers. At the lowest pressures (30 and 50 GPa), the bottom and top layer protons do not mix: although proton hopping (either of a few beads or of the whole ring-polymer in the PIMD framework) from one layer to the other does occur, the protons always return to their original layer by a second hop in reverse. Protons initially situated on one layer thus remain on that layer for the whole duration of the simulation with, from time to time, a short exploration of the other layer. However, at 70 GPa, lower and upper layers are not distinguishable as reverse hopping does not always follow: protons do not belong to one particular layer anymore but rather all protons should be considered as forming a single quasi-2D plane within which efficient diffusion can occur, as we further analyze in the following. Interestingly at the highest probed pressure, 90 GPa, the overall distribution becomes narrower, but the lower and upper protons can again be distinguished: as for the lowest two pressures, they tend to stay in their initial layer with only occasional hops, despite a large overlap between the partial distributions. This behavior will be rationalized in the following, examining the different pressure trends of the two hopping mechanisms that are necessary for long-range proton diffusion.

4/12
In-plane reorientation First we discuss the reorientation mechanism. As already described, within the P3 structure, protons hop in-plane between the 6i sites. This motion can be efficiently described by the azimuthal angle φ as shown in the sketch of figure 3. From the probability distribution of the latter P r (φ ), we extracted the Gibbs free-energy profile G = −k B T log P r (φ ), which includes both thermal and quantum effects.
As shown in Figure 3, the proton free-energy profile along this coordinate has a three-fold symmetry with equivalent barrier heights between the three wells, consistently with the P3 configuration. The barrier width of this free energy profile is grossly proportional to 2π Upon compression, we observe that the free energy barrier height increases, from 20 meV at 30 GPa up to 100 meV at 90GPa, revealing a pressure induced confinement of the proton along this coordinate. This is in line with the increasing of the average polar angle θ with pressure that splits the 6i sites apart. Along with the compression of the layers along the c axis, the protons from the two layers try to minimize their mutual repulsion by increasing θ . As a consequence, the reorientational dynamical disorder, thermally activated at low pressure, tends to slow down.We point out that classical simulations, not including NQE, yield almost identical P r (φ ) distributions, meaning that the quantum behavior is in this case moderate within the pressure range that we explored. The effect of pressure contrasts with that of temperature, which tends to allow the proton to explore equivalently the three wells 33 by thermal activation.

Out-of-plane dissociation
As the hydrogen planes become closer upon compression, the protons in 6i sites adopt quasi-linear O-H-O configurations, where the O anions belong to distinct Mg(OH) 2 stacks. This configuration is thus prone to the formation of inter-layer hydrogen bonds. It has been shown 13 that only weak hydrogen bonds could be present in brucite. However, as we discuss below, a double-well potential is found along the O-O direction at moderate pressure, which suggests the onset of hydrogen bonds and occurs in parallel with the pressure-induced creation of quasi-2D proton layers in the structure. All along those transformations, NQE play an important role.
In order to investigate the proton effective potential, we adopt as the order parameter 18 χ which is the difference between the distances that separate the hydrogen atoms from their nearest and second nearest neighbor oxygen atoms projected on the O-O direction (see sketch in figure 4). bringing the two proton equilibrium positions closer and reducing |χ 0 |. At high pressures, a proton can therefore hop from one oxygen atom to another through quantum tunneling, zero-point energy or thermal activation, which constitutes the so-called "dissociation" process 33 . The difference between classical and quantum simulations is significant as the classical barrier is ∼ 3k B T higher than its quantum counterpart at 70 GPa. Therefore, while the proton reorientation mechanism is thermally activated, the dissociation process is mainly quantum driven.

Proton diffusion sweet spot
The evolution of the free energy barrier heights upon compression is reported in figure 5. In the case of brucite Mg(OH) 2 , the barrier height for dissociation decreases from ∆G d ∼ 0.11 eV at 30 GPa to ∆G d ∼ 0.01 eV at 90 GPa. In contrast, the reorientation barrier increases from ∆G r ∼ 0.03 eV to ∆G r ∼ 0.1 eV within the same pressure range. The two curves cross at ∼ 70 GPa where hopping rates for the two processes have the same order of magnitude, while reorientation dominates at lower pressures and dissociation prevails at higher pressures. We thus suggest that the highest proton mobility occurs around P = 70 GPa, which represents the sweet spot for proton diffusion. We provide in table 2 a rough estimate of both dissociation and reorientation reaction rates κ d and κ r , through the Eyring-Polanyi equation 49,50 : While in brucite the typical reorientation and dissociation times are reachable within our simulation duration, in portlandite (discussed in greater detail below), because of the larger barriers with respect to brucite (see Fig.5), proton dissociation characteristic time in Ca(OH) 2 is on the ns scale or even larger for all pressures before the amorphization transition. We note in passing that κ −1 , as evaluated for classical protons within the Eyring-Polanyi equation, is several order of magnitudes greater.
Brucite Mg(OH) 2 Portlandite Ca(OH) 2 P (GPa) 30   Free energy barriers for brucite and portlandite, as computed from the probability distributions for reorientation (∆G r ) and dissociation (∆G d ). As discussed in the text, the barriers for proton reorientation and diffusion show opposite trends with increasing pressure, with a non-linear increase of the reorientation barrier in Mg(OH) 2 . The competition between the two trends generates a pressure sweet spot in brucite between 60 and 70 GPa, where both reorientation and dissociation barriers can be overcome; in this pressure range, proton diffusion is enhanced, which requires both mechanisms to occur. Such a sweet spot cannot be reached in portlandite due to the transition taking place at 6 GPa.

In-plane proton distribution.
The barrier height analysis above is perfectly consistent with the probability distribution of the proton positions in the (a, b) plane as shown in figure 6. For P = 30 GPa, the proton distribution shows three broad peaks next to each oxygen atom, thus revealing reorientation processes between the 6i sites. As pressure is increased to 50 GPa, the peaks become narrower and centered at a greater distance from the oxygen sites. This is due to the hindering of the reorientations and the progressive slanting of the O-H bonds towards the (a, b) plane. In addition, the proton density midway between the oxygen sites increases as dissociation become easier. At P = 70 GPa, we observe evidence of a long-range proton diffusion process as the hydrogen nuclei spread all over the simulation box. The onset of dissociation, while reorientation processes still occur, allows the protons to migrate beyond their first neighbors. The range of this in-plane motion (over several Å) can be compared with the width of the out-of-plane c axis distribution as in figure 2 (approximately 0.2Å), clearly indicating the two-dimensional character of this process. Finally, at the highest pressure, P = 90 GPa, although the dissociation probability keeps increasing, the reorientation processes are locked in. The computed proton distribution is not symmetric among the 6i sites anymore and long-range proton diffusion is hindered. This confirms the specificity of the 60-70 GPa pressure range as a sweet spot for proton diffusion.
Due to the relative shortness of the PIMD runs and the fact that a Langevin thermostat with a friction coefficient of 10ps −1 is attached to the centroid motion to achieve efficient sampling (see Methods section), the PIMD simulations cannot be used to obtain a direct estimate of the proton diffusion coefficient. Nonetheless, the trends displayed on Fig. 6, with the onset of efficient proton diffusion at the sweet spot pressure of 70GPa are qualitatively correct and consistent with the previous discussion. Furthermore, the PILE-L thermostatting scheme 36 that we use in this work reduces, in the limit of a vanishing centroid friction, to the Thermostatted Ring-Polymer Molecular Dynamics (TRPMD) algorithm 51 for computing dynamical observables. We expect centroid thermostatting to hinder proton diffusion in our PIMD dynamics, which explains why long-range proton diffusion is not directly observed at 50GPa, even though the estimated dissociation and reorientation inverse rates are relatively short compared to the PIMD trajectory length. Note that centroid thermostatting does not affect the PIMD estimates for static properties such as the free energy barrier ∆G, and therefore it has no impact on the rates presented in Table 2. Proton diffusion should indeed occur already at pressures lower than 70GPa and it might be at the root of a de-stabilization of the system easing the suggested 42 phase transition, or could be involved in the water transfer in the earth mantle.

Comparison with portlandite
Finally, we close our discussion on proton diffusion in brucite-like minerals by a comparison with portlandite Ca(OH) 2 which has the brucite structure for pressures up to approximately 6 GPa. The same analysis as for brucite was done systematically for  Figure 6. Distribution of the proton positions in the (x, y) plane at 30, 50, 70 and 90 GPa for brucite Mg(OH) 2 (upper panels) and portlandite Ca(OH) 2 (lower panels). The circles represent the projection of the oxygen sites on the (x, y) plane (light red: bottom layer, light blue: top layer). Periodic boundary conditions are not used in computing this distribution in order to visualize the displacement of the protons: in brucite, the protons that were initially located within the simulation box move out of its boundaries at 70 GPa, which is an indication of diffusion occurring even within the 30 ps duration of our simulations. In portlandite, this might happen at pressures above 15 GPa, thus beyond the structural phase transition pressure reported around 6 GPa from refs. 40, 41 . portlandite. In figure 5 we present the evolution of free energy barriers of the proton reorientation and dissociation mechanisms. We observe that in portlandite the reorientational barrier ∆G r at 10 GPa is comparable to that in brucite at 50 GPa. However, the pressure effect on the latter barrier is more important in portlandite as shown by the larger pressure slope of 2.8 meV/GPa while its counterpart in brucite is about 1.6 meV/GPa. Ca ++ cations are larger than Mg ++ ones, so that the computed in-plane O-O distance at P = 10 GPa is 3.45 Å in portlandite and 3.06 Å in brucite; in contrast, the distance between the oxygen anions on distinct stacks is almost the same in the two crystals. This implies larger polar angles θ in portlandite than in brucite, which efficiently hinder the reorientation mechanism, even at relatively low pressures. As far as the dissociation barriers ∆G d are concerned, in portlandite they are greater than in brucite but decrease much faster with pressure. The estimated pressure slope ∂ ∆G d ∂ P is 14 meV/GPa in portlandite, as compared to 2meV/GPa in brucite. This derives from the larger compressibility of portlandite with respect to brucite, as demonstrated in a recent work 11 . It has to be noticed that the large value of the dissociation barrier in portlandite would require very long runs of path integral ab initio molecular dynamics, much beyond the scope of this work, to record a significant number of events. Nevertheless, we detected a few events as some of the replicas of the PIMD simulations did occasionally reach the barrier top. A precise evaluation of the barrier heights would have nevertheless required the use of accelerated sampling techniques, such as in ref. 52.
Finally, the crossing point of the dissociation and orientational barriers in portlandite should occur beyond 20 GPa, with diffusion barriers comparable to that of brucite at 70 GPa. However, a transition towards another phase is reported at 6GPa 40,41 and our own simulations reveal instability of the system at 20 GPa. Therefore, as shown in figure 6, no diffusion was observed for portlandite within the time scale of our simulations. Indeed, the reaction rate estimates, given in table 2, yield much longer times than in brucite.
The values reported in Table 2 for portlandite agree remarkably well with the results of Dupuis and coworkers 33 as regards the reorientation rate at 10 GPa and its increase with decreasing pressure. Our results for the dissociation rate, however, are lower than that of ref. 33 by as much as 3 orders of magnitude under comparable thermodynamic conditions. In ref. 33, different methodologies were used for the computation of the two rates: whereas path integral methods were employed for 8/12 the reorientation process, metadynamics with classical nuclei was used for dissociation. In this work on the other hand, we compute the two rates within a fully consistent framework: the same PIMD trajectories are employed to estimate the free energy barriers associated with the two processes, which are then inserted into the Eyring-Polanyi approximation. We do not resort to any type of accelerated sampling scheme to overcome the barriers. Note that for the lowest pressure of 5 GPa, the ring-polymer beads very rarely explore the top of the dissociation barrier, therefore the uncertainty on ∆G d is large (we estimate ∆G d = ± 20 meV). This accuracy is nonetheless sufficient to exclude, within our approach, such a fast dissociation as reported in ref. 33.
It can be noted, that among the other systems sharing the same structure as brucite, our preliminary analysis of theophrastite Ni(OH) 2 indicates that this system should also present a sweet spot for proton diffusion at approximately the same pressure as Mg(OH) 2 , due to comparable ionic radii between Mg ++ and Ni ++ cations.

Conclusion
To summarize, we analyzed the proton diffusion mechanism in both brucite (Mg(OH) 2 ) and portlandite (Ca(OH) 2 ) under pressure, taking into account nuclear quantum effects by path-integral based ab initio molecular dynamics. Proton diffusion in those crystals involves two stages to occur: a reorientation motion within the (a,b) plane, and a proton dissociation between two oxygen atoms on opposite layers. Firstly, we have seen that the reorientation mechanism is thermally activated and that the pressure tends to localize the proton in a certain orientation, making the reorientation motion less likely. Secondly, in contrast with the reorientation, we showed that the dissociation mechanism was quantum driven and was made easier by increasing pressure through the formation of a quasi two-dimensional hydrogen layer.
These two antagonistic effects give rise to a pressure sweet spot for proton diffusion within 60 and 70 GPa in brucite. However, proton diffusion could also occur at much lower pressure, although it is less probable, and could be at the root of a destabilization of the structure, as suggested by the theoretical predictions of a phase transition 42 at 20 GPa or decomposition into MgO and H 2 O at the same pressure 6 . Beyond this pressure threshold the reorientation becomes a bottleneck for proton diffusion, while dissociation is the rate-limiting step at lower pressure.
Finally, by systematic comparison with portlandite, we demonstrate the specificity of brucite for proton diffusion. Indeed, the proton diffusion sweet spot in portlandite would occur at pressures well beyond its transition toward an another phase.

Methods
Molecular dynamics (MD) simulations were carried out at room temperature and fixed volume either via a classical Langevin equation or within the Path Integral (PI) framework to take into account nuclear quantum effects while the electronic structure and atomic forces were described within the Density Functional Theory (DFT). We compute the electron density and atomic forces within the Perdew-Wang Generalized Gradient approximation to the exchange and correlation density functional 53 , using the Quantum Espresso package 54 in combination with the i-PI 55 interface for the path-integral simulations. We employed ultra-soft pseudo-potentials with a plane wave expansion cutoff of E cut = 50 Ry for the Kohn-Sham states and 8 times as large for the charge and the potential, ensuring total energy convergence. Both brucite Mg(OH) 2 and portlandite Ca(OH) 2 were simulated by using a hexagonal ( √ 3 × √ 3 × 1) supercell and a (2 × 2 × 2) k-point sampling centred at 2π a ( 1 2 , 1 2 , a 2c ). The number of beads in the PIMD simulations was set to 24 and checked to provide convergence of kinetic and potential energies. The PIMD simulations were performed using the efficient PILE-L thermostatting scheme 36 with a centroid friction coefficient of 10 ps −1 . Lattice parameters were obtained through systematic volume relaxation of the system ensuring isotropic stress tensors for each pressure. The optimized equilibrium lattice parameters a and c were remarkably similar between classical and path-integral MD simulations for the pressures here considered: the differences between the classical and quantum results for a and c amounted to few hundredths of Å. Finally, the typical duration time of the simulations was 30 ps.