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

Sofiane Schaack, Philippe Depondt, Simon Huppert, Fabio Finocchi

► To cite this version:

Sofiane Schaack, Philippe Depondt, Simon Huppert, Fabio Finocchi. Quantum driven proton diffusion in brucite-like minerals under high pressure. Scientific Reports, 2020, 10 (1), 10.1038/s41598-020-64813-8. http://dx.doi.org/10.1038/s41598-020-64813-8

HAL Id: hal-02935811

https://hal.sorbonne-universite.fr/hal-02935811

Submitted on 10 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Sofiane Schaack, Philippe Depondt,*, Simon Huppert & Fabio Finocchi

Transport of hydrogen in hydrous minerals under high pressure is a key step for the water cycle within the Earth interior. Brucite Mg(OH)₂ is one of the simplest minerals containing hydroxyl groups and is believed to decompose under the geological condition of the deep Earth's mantle. In the present study, we investigate the proton diffusion in brucite under high pressure, which results from a complex interplay between two processes: the O-H reorientations motion around the 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 Mg(OH)₂, proton diffusion reaches a maximum for pressures close to 70GPa, while the structurally similar portlandite Ca(OH)₂ never shows proton diffusion within the pressure range and time scale that we explored. We analyze the different behavior of brucite and portlandite, which might constitute two prototypes for other minerals with same structure.

Hydroxide minerals play an important role in several problems in geology, surface science or for industrial applications. Among them, brucite Mg(OH)₂ can be formed at the interface between periclase MgO and water at ambient conditions¹⁻³. The trigonal brucite structure consists of alternating layers along the c axis that terminate with hydroxyls (Fig. 1). This structure is common to other hydroxides of divalent metals, such as Ca(OH)₂, Ni(OH)₂ and Cd(OH)₂. Portlandite Ca(OH)₂ 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)₂ can also act as a water vector in subduction zones, through complex processes that take place within the Earth interior⁴,⁵. Therefore, the behavior of brucite and brucite-like minerals at very high pressure has been widely investigated. X-ray diffraction of Mg(OH)₂ up to 78 GPa showed that the c/a ratio decreases steadily from ambient pressure up to about 25 GPa and then stays almost constant⁶. 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 parallel to the c axis at ambient conditions slant in three equivalent positions as the inter-layer distance shrinks even under moderate pressure⁷ or when decreasing temperature⁸. Besides reducing the global symmetry from P3$m_1$ 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⁹, which eventually reinforce under further compression and modify the compressibility along c¹⁰. 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¹¹,¹², has also been invoked as the reason for the pressure-induced hydrogen sublattice amorphization in brucites⁷,¹³. 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, that is, all the properties that go beyond a purely classical description of ion dynamics¹⁴, such as zero-point

Institut des NanoSciences de Paris (INSP), Sorbonne Université, CNRS-UMR 7588, 75005, Paris, France. e-mail: depondt@insp.jussieu.fr
<figure><img src="image_3.png" /><figcaption>**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)₂ units. The hydrogen atoms in the P₃m₁ symmetry structure (light grey) are in the 2d Wyckoff sites, while the hydrogen atoms in the P₃ 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.</figcaption></figure>

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¹⁵⁻¹⁷.

Despite the complexity and computational cost of accounting for the quantum nature of light nuclei in simulations, efficient methods have arisen¹⁸⁻²¹, and, stimulated by the increasing interest in often paradoxical nuclear quantum effects (NQE)¹⁴, a new field is rapidly growing with potential applications in an eclectic collection of issues including protonic conduction²², hydrogen in biological matter²³, water circulation in Earth's mantle²⁴,²⁵ or hydrogen storage²⁶,²⁷.

In this work, we address the proton diffusion process occurring in brucite Mg(OH)₂ that appears to be a specific case as compared with portlandite Ca(OH)₂, which was the object of a recent theoretical study²⁸, and is also discussed here. The reasons for the different behaviors of these two compounds are elucidated. We employ classical and Path-Integral Molecular Dynamics²⁹⁻³¹, (PIMD) to disentangle thermal from nuclear quantum effects in a consistent framework. Among the various proton diffusion mechanisms that were considered in the past, we confirm that the *relay mechanism*²⁸,³², is at play in brucite under high pressure. This two-step process involves proton hops from one oxygen atom to another and rotations of OH groups. We also show in the following that proton hopping is quantum driven while OH rotation is mostly thermally activated.

## Structural Changes Under Pressure

The brucite minerals are layered structures composed by stacks of metallic ions, oxygen, and hydrogen atoms in a CaI₂-type structure¹¹,¹²,²⁸. The metallic element (Mg, Ca, Ni...) has an impact on several physical parameters, in particular, the lattice parameters and the compressibility of the system are different in the two systems we studied: brucite (Mg(OH)₂) and portlandite (Ca(OH)₂)¹⁰. The metal layers terminate with hydroxyl groups on both sides (Fig. 1, left). At ambient pressure and temperature, the brucite minerals belong to the P₃m₁ space group with hydrogen atoms located on the threefold axis above or below oxygen atoms (2d Wyckoff sites) with full occupancy as shown in the left-hand side of Fig. 1.

However, in accordance with previous experimental and simulation results⁷,¹¹,¹², as pressure increases, due to the growing repulsive interaction between protons on opposite layers, the protons do not remain above the corresponding oxygen atom. Brucite minerals adopt a P₃ configuration, where the H nuclei are in the 6i Wyckoff sites with a 1/3 occupancy factor (Fig. 1, left). This proton rearrangement has two noticeable consequences: first, a frustration of the proton orientation upon compression¹¹,³³,³⁴ second, the protons in 6i sites could form hydrogen bonds between distinct layers⁷. This point, which has been raised by several authors, is discussed in the following.

The stability of both brucite and portlandite upon compression has been investigated in the past and an amorphization of portlandite was found through infra-red absorption to occur between 10 and 15 GPa³³. In brucite, however, X-ray experiments⁶ up to 78 GPa at 600 K showed a smooth variation of the diffraction patterns; in the same paper, the authors deduced an equation of state according to which brucite would decompose into periclase
and water around 27 GPa at ambient temperature⁶. These apparently contradictory results suggest a significant kinetic decomposition barrier. More recently, a phase transition toward a tetragonal structure was suggested in brucite above 20 GPa by ab initio simulations³⁵; nonetheless, up to now, no such transition was observed experimentally to our knowledge⁶,³⁴. The pressure increase shrinks the distance between the MgO layers, due to the important compressibility¹⁰ along the c axis, thus allowing 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 to decompose. Moreover, the formation of inter-layer hydrogen bonds could allow the protons to hop between the two oxygen atoms on the facing Mg(OH)₂ layers.

**Proton diffusion mechanism.** As in other hydrous minerals²⁴,²⁵, hydrates³⁶ and ices¹⁵,³⁷, the proton can hop along hydrogen bonds. Following Dupuis et al.²⁸, 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 P 3̄ 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³²,³⁸.

In P3̄ Mg(OH)₂, 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²⁸ 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)₂ and Ca(OH)₂ 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.

**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 $\phi$ as shown in the sketch of Fig. 3. From the probability distribution of the latter $\mathcal{P}_t(\phi)$, we extracted the Gibbs free-energy profile $G = -k_B T \log \mathcal{P}_t(\phi)$, which includes both thermal and quantum effects.

As shown in Fig. 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 $\frac{2\pi}{3} d_{O-H} \cos\theta$ ($\theta$ being the zenith polar angle and $d_{O-H}$ the covalent O-H bond length), that is, how far the O-H bond slants away from the c axis.

Upon compression, we observe that the free energy barrier height increases, from 20 meV at 30 GPa up to 100 meV at 90 GPa, revealing a pressure induced confinement of the proton along this coordinate. This is in line
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Probability distribution of the proton position along *c*. The green lines relate to the bottom layer of hydrogen nuclei (Fig. 1), the blue lines to the top layer, while the purple lines are the sum of both.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Free energy profile along the polar angle $\phi$. The proximity between classical and PIMD profiles at 30 GPa shows that the reorientation mechanism is dominated by thermal effects. The thermal energies corresponding to $T = 300$ K and higher temperatures are also shown as long-dotted lines for comparison.</figcaption></figure>

with the increasing of the average polar angle $\theta$ 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 $\theta$. 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_\phi(\phi)$ 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²⁸ 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)₂ stacks. This configuration is thus prone to the formation of inter-layer hydrogen bonds. It has been shown¹² that only weak
<figure><img src="image_6.png" /><figcaption>Figure 4. Proton free energy profile along the order parameter χ as defined in the text and sketched on the right (O in red, proton in gray).</figcaption></figure>

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¹⁷ χ 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 Fig. 4).

$$ \chi = |\vec{O_2H} \cdot \vec{u_{OO}}| - |\vec{O_1H} \cdot \vec{u_{OO}}| \quad (1) $$

with $\vec{u_{OO}}$ the unitary vector in the O₁-O₂ direction.

The free-energy profiles along this coordinate (see Fig. 4) show that: (i) the effective potential is a double well with minima at $\pm \chi_0$, as for other hydrogen-bond forming crystals¹⁷,²⁴; (ii) the barrier height decreases upon compression, from 4 k_BT at P = 30 GPa to 0.5 k_BT at P = 90 GPa. This occurs while the O-O distance shrinks with pressure along the c axis, bringing the two proton equilibrium positions closer and reducing $|\chi_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²⁸. The difference between classical and quantum simulations is significant as the classical barrier is ~3k_BT 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 Fig. 5. In the case of brucite Mg(OH)₂, the barrier height for dissociation decreases from $\Delta G_d \sim 0.11$ eV at 30 GPa to $\Delta G_d \sim 0.01$ eV at 90 GPa. In contrast, the reorientation barrier increases from $\Delta G_r \sim 0.03$ eV to $\Delta G_r \sim 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 1 a rough estimate of both dissociation and reorientation reaction rates $\kappa_d$ and $\kappa_r$, through the Eyring-Polanyi equation³⁹,⁴⁰:

$$ \kappa_{d,r} = \frac{k_B T}{h} e^{-\frac{\Delta G_{d,r}}{k_B T}} \quad (2) $$

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)₂ is on the ns scale or even larger for all pressures before the amorphization transition. We note in passing that $\kappa^{-1}$, as evaluated for classical protons within the Eyring-Polanyi equation, is several order of magnitudes greater.

**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 Fig. 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
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Free energy barriers for brucite and portlandite, as computed from the probability distributions for reorientation (ΔG<sub>r</sub>) and dissociation (ΔG<sub>d</sub>). 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)<sub>2</sub>. 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 amorphization transition taking place between 10 and 15 GPa.</figcaption></figure>

<table><caption>**Table 1.** Eyring-Polanyi inverse reaction rate for the dissociation $κ_d^{-1}$ and reorientation $κ_r^{-1}$ mechanisms, computed for brucite (left) and portlandite (right). The temperature is T = 300 K. The transition pressure for amorphization of portlandite is reported from refs. 13,42.</caption><thead><tr><th></th><th></th><th colspan="4">Brucite Mg(OH)<sub>2</sub></th><th colspan="2">Portlandite Ca(OH)<sub>2</sub></th><th></th><th></th></tr><tr><th>P</th><th>(GPa)</th><th>30</th><th>50</th><th>70</th><th>90</th><th>5</th><th>10</th><th></th><th>15</th></tr></thead><tbody><tr><td>κ<sub>r</sub><sup>-1</sup></td><td>(ps)</td><td>0.43</td><td>0.65</td><td>1.47</td><td>12.3</td><td>0.32</td><td>0.54</td><td rowspan="2">Transition</td><td>0.95</td></tr><tr><td>κ<sub>d</sub><sup>-1</sup></td><td>(ps)</td><td>13.2</td><td>3.07</td><td>0.83</td><td>0.24</td><td>8.63 × 10<sup>3</sup></td><td>438</td><td>44.3</td></tr></tbody></table>

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 Fig. 2 (approximatively 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 10 ps⁻¹ 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 70 GPa are qualitatively correct and consistent with the previous discussion. Furthermore, the PILE-L thermostatting scheme³¹ that we use in this work reduces, in the limit of a vanishing centroid friction, to the Thermostatted Ring-Polymer Molecular Dynamics (TRPMD) algorithm⁴¹ 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 50 GPa, 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 1. Proton diffusion should indeed occur already at pressures lower than 70 GPa and it might be at the root of a de-stabilization of the system easing the suggested³⁵ 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)<sub>2</sub> which has the brucite structure for pressures up to approximately 15 GPa. The same analysis as for brucite was done systematically for portlandite. In Fig. 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<sub>r</sub> 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<sup>++</sup> cations are larger than Mg<sup>++</sup> ones, so that the computed
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Distribution of the proton positions in the (x, y) plane at 30, 50, 70 and 90 GPa for brucite Mg(OH)₂ (upper panels) and portlandite Ca(OH)₂ (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 amorphization transition pressure that is evidenced by the red vertical line.</figcaption></figure>

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 dθ/dG_d/dp is 14 meV/GPa in portlandite, as compared to 2 meV/GPa in brucite. This derives from the larger compressibility of portlandite with respect to brucite, as demonstrated in a recent work¹⁰. 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. ⁴².

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 an amorphous phase is reported between 10 and 15 GPa¹³,⁴³, and our own simulations reveal instability of the system at 20 GPa. Therefore, as shown in Fig. 6, no diffusion was observed for portlandite within the time scale of our simulations. Indeed, the reaction rate estimates, given in Table 1, yield much longer times than in brucite.

The values reported in Table 1 for portlandite agree remarkably well with the results of Dupuis and coworkers²⁸ 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. ²⁸ by as much as 3 orders of magnitude under comparable thermodynamic conditions. In ref. ²⁸, different methodologies were used for the computation of the two rates: whereas path integral methods were employed for 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. ²⁸.
It can be noted, that among the other systems sharing the same structure as brucite, our preliminary analysis of theophrastite Ni(OH)₂ indicates that this system should also present a sweet spot for proton diffusion at approximately the same pressure as Mg(OH)₂, due to comparable ionic radii between Mg⁺⁺ and Ni⁺⁺ cations.

## Conclusion

To summarize, we analyzed the proton diffusion mechanism in both brucite (Mg(OH)₂) and portlandite (Ca(OH)₂) 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³⁵ at 20 GPa or decomposition into MgO and H₂O at the same pressure⁶. 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 amorphous 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¹⁴, using the Quantum Espresso package⁴⁵ in combination with the i-PI⁴⁶ 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)₂ and portlandite Ca(OH)₂ were simulated by using a hexagonal $(\sqrt{3} \times \sqrt{3} \times 1)$ supercell and a $(2 \times 2 \times 2)$ k-point sampling centred at $(\frac{2\pi}{a}(\frac{1}{2}, \frac{1}{2}, \frac{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³¹ with a centroid friction coefficient of 10 ps⁻¹. 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.

Received: 14 November 2019; Accepted: 3 February 2020;

Published online: 15 May 2020

## References

1. Chemical weathering rates of silicate mineral. In White, A. & Brantley, S. (eds.) Reviews in Mineralogy, vol. 31 (Mineralogical Society of America, 1995).

2. Jordan, G., Higgins, S. & Egglestone, C. Dissolution of the periclase (001) surface: a scanning force microscope study. *Am. Miner.* **84** (1999).

3. Sharma, R., McKelvy, M. J., Béarat, H., Chizmeshya, A. V. & Carpenter, R. In-situ nanoscale observations of the Mg(OH)₂ dehydroxylation and rehydroxylation mechanisms. *Philosophical Magazine* **84**, 2711–2729 (2004).

4. Bodnar, R. et al. Whole earth geohydrologic cycle, from the clouds to the core: The distribution of water in the dynamic earth system. In Bickford, M. (ed.) *The Web of Geological Sciences: Advances, Impacts, and Interactions*, vol. Special Papers (Geological Society of America, 2013).

5. Peacock, S. & Hyndman, R. Hydrous minerals in the mantle wedge and the maximum depth of subduction thrust earthquakes. *Geophys. Res. Lett.* **26** (1999).

6. Fei, Y. & Mao, H.-K. Static compression of Mg(OH)₂ to 78 GPa at high temperature and constraints on the equation of state of fluid H₂O. *Journal of Geophysical Research: Solid Earth* **98**, 11875–11884 (1993).

7. Catti, M., Ferraris, G., Hull, S. & Pavese, A. Static compression and h disorder in brucite, Mg(OH)₂, to 11 gpa: a powder neutron diffraction study. *Physics and Chemistry of Minerals* **22**, 200–206 (1995).

8. Chakoumakos, B., Loong, C.-K. & Schultz, A. Low-temperature structure and dynamics of brucite. *The Journal of Physical Chemistry B* **101**, 9458–9462 (1997).

9. Parise, J. B., Leinenweber, K., Weidner, D. J., Tan, K. & Dreele, R. B. V. Pressure-induced h bonding: Neutron diffraction study of brucite, Mg(OH)₂, to 9.3 GPa. *American Mineralogist* **79**, 193–196 (1994).

10. Ulian, G. & Valdrè, G. Equation of state and second-order elastic constants of portlandite Ca(OH)₂ and brucite Mg(OH)₂. *Physics and Chemistry of Minerals* **46**, 101–117 (2019).

11. Raugei, S., Silvestrelli, P. L. & Parrinello, M. Pressure-induced frustration and disorder in Mg(OH)₂ and Ca(OH)₂. *Physical Review Letters* **83**, 2222 (1999).

12. Mookherjee, M. & Stixrude, L. High-pressure proton disorder in brucite. *American Mineralogist* **91**, 127–134 (2006).

13. Nagai, T., Hattori, T. & Yamanaka, T. Compression mechanism of brucite: An investigation by structural refinement under pressure. *American Mineralogist* **85**, 760–764 (2000).

14. Markland, T. E. & Ceriotti, M. Nuclear quantum effects enter the mainstream. *Nature Rev. Chem.* **2** (2018).

15. Benoit, M., Marx, D. & Parrinello, M. Tunnelling and zero-point motion in high-pressure ice. *Nature* **258–261** (1998).
16. Dammak, H., Antoshchenkova, E., Hayoun, M. & Finocchi, F. Isotope effects in lithium hydride and lithium deuteride crystals by molecular dynamics simulations. *Journal of Physics: Condensed Matter* **24**, 435402 (2012).

17. Bronstein, Y., Depondt, P., Finocchi, F. & Saitta, A. M. Quantum-driven phase transition in ice described via an efficient langevin approach. *Phys. Rev. B* **89**, 214101 (2014).

18. Berne, B. & Thumalai, R. On the simulation of quantum systems: path integral methods. *Annual Review of Physical Chemistry* **37** (1986).

19. Marx, D. & Parrinello, M. *Ab initio* path integral molecular dynamics: Basic ideas. *J. Chem. Phys.* **104**, 4077 (1996).

20. Tuckerman, M. E., Marx, D., Klein, M. L. & Parrinello, M. Efficient and general algorithms for path integral Car - Parrinello molecular dynamics. *J. Chem. Phys.* **104**, 5579 (1996).

21. Dammak, H., Chalopin, Y., Laroche, M., Hayoun, M. & Greffet, J.-J. Quantum thermal bath for molecular dynamics simulation. *Phys. Rev. Lett.* **103**, 190601 (2009).

22. Liang, K. C. & Nowick, A. S. High-temperature protonic conduction in mixed perovskite ceramics. *Solid State Ionics* **61**, 77–81 (1993).

23. Rossi, M., Ceriotti, M. & Manolopoulos, D. E. Nuclear quantum effects in H⁺ and OH⁻ diffusion along confined water wires. *The Journal of Physical Chemistry Letters* **7**, 3001–3007 (2016).

24. Bronstein, Y., Depondt, P. & Finocchi, F. Thermal and nuclear quantum effects in the hydrogen bond dynamical symmetrization phase transition of δ-AlOOH. *European Journal of Mineralogy* **29**, 1 (2017).

25. Sano-Furukawa, A. et al. Direct observation of symmetrization of hydrogen bond in δ-AlOOH under mantle conditions using neutron diffraction. *Scientific Reports* **8**, 15520 (2018).

26. Zhang, C. & Michaelides, A. Quantum nuclear effects on the location of hydrogen above and below the palladium (100) surface. *Surface Science* **605**, 689–694 (2011).

27. Wahiduzzaman, M., Walther, C. F. J. & Heine, T. Hydrogen adsorption in metal-organic frameworks: The role of nuclear quantum effects. *The Journal of Chemical Physics* **141**, 064708 (2014).

28. Dupuis, R., Dolado, J. S., Benoit, M., Surga, J. & Ayuela, A. Quantum nuclear dynamics of protons within layered hydroxides at high pressure. *Scientific reports* **7**, 4842 (2017).

29. Feynman, R. P., Hibbs, A. R. &Styer, D. F. *Quantum mechanics and path integrals* (Courier Corporation, 2010).

30. Chandler, D. & Wolynes, P. G. Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids. *J. Chem. Phys.* **74**, 4078–4095 (1981).

31. Ceriotti, M., Parrinello, M., Markland, T. E. & Manolopoulos, D. E. Efficient stochastic thermostatting of path integral molecular dynamics. *J. Chem. Phys.* **133**, 124104, https://doi.org/10.1063/1.3489925 (2010).

32. Yaroslavstev, A. B. Proton conduction of inorganic hydrates. *Russian Chemical Reviews* **63**, 429–435 (1994).

33. Kruger, M., Williams, Q. & Jeanloz, R. Vibrational spectra of Mg(OH)₂ and Ca(OH)₂ 2 under pressure. *The Journal of Chemical Physics* **91**, 5910–5915 (1989).

34. Duffy, T. S., Meade, C., Fei, Y., Mao, H.-K. & Hemley, R. J. High-pressure phase transition in bruciteMg(OH)₂. *American Mineralogist* **80**, 222–230 (1995).

35. Hermann, A. & Mookherjee, M. High-pressure phase of brucite stable at earths mantle transition zone and lower mantle conditions. *Proceedings of the National Academy of Sciences* **113**, 13971–13976 (2016).

36. Schaack, S. et al. Observation of methane filled hexagonal ice stable up to 150 GPa. *Proc. Natl. Acad. Sci. USA* **116**, 16204 (2019).

37. Bronstein, Y. et al. Quantum versus classical protons in pure and salty ice under pressure. *Phys. Rev. B* **93**, 024104 (2016).

38. de Grotthuss, C. T. Sur la décomposition de l'eau et des corps qu'elle tient en dissolution À l'aide de l'électricité galvanique. *Annales de Chimie* **58**, 54–73 (1806).

39. Eyring, H. The activated complex in chemical reactions. *J. Chem. Phys.* **3** (1935).

40. Evans, M. & Polanyi, M. Some applications of the transition state method to the calculation of reaction velocities, especially in solution. *Trans. Faraday Soc.* **31** (1935).

41. Rossi, M., Ceriotti, M. & Manolopoulos, D. E. How to remove the spurious resonances from ring polymer molecular dynamics. *J. Chem. Phys.* **140**, 4116 (2014).

42. Ivanov, S., Grant, I. &Marx, D. Quantum free energy landscapes from ab initio path integral metadynamics: Double proton transfer in the formic acid dimer is concerted but not correlated. *J. Chem. Phys.* **143** (2015).

43. Meade, C. & Jeanloz, R. Static compression of Ca(OH)₂ at room temperature: observations of amorphization and equation of state measurements to 10.7 GPa. *Geophysical Research Letters* **17**, 1157–1160 (1990).

44. Burke, K., Perdew, J. &Wang, Y. Derivation of a generalized gradient approximation: the pw91 density functional. In Dobson, J., Vignale, G. &Das, M. (eds.) *Electronic Density Functional Theory: recent progress and new directions* (Springer Science, 1996).

45. Giannozzi, P. et al. Quantum espresso: a modular and open-source software project for quantum simulations of materials. *J. Phys.: Cond. Matt* **21**, 395502 (2009).

46. Ceriotti, M., More, J. &Manolopoulos, D. E. i-pi: A python interface for ab initio path integral molecular dynamics simulations. *Computer Phys. Comm.* **189** (2014).

## Acknowledgements

This work was granted access to the HPC resources of CINES under the allocation A0030906719 made by GENCI.

## Author contributions

S.S., P.D. and F.F. designed research; S.S. carried out the simulations; S.S., P.D., and S.H. analyzed the results; S.S. and P.D. wrote the manuscript that was jointly revised by all of the authors.

## Competing interests

The authors declare no competing interests.

## Additional information

Correspondence and requests for materials should be addressed to P.D.

Reprints and permissions information is available at www.nature.com/reprints.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or
format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Cre-
ative Commons license, and indicate if changes were made. The images or other third party material in this
article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the
material. If material is not included in the article's Creative Commons license and your intended use is not per-
mitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the
copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020