Ab initio investigation of H-bond disordering in δ -AlOOH

δ -AlOOH ( δ ) is a high-pressure hydrous phase that participates in the deep geological water cycle. At 0 GPa, δ has asymmetric hydrogen bonds (H-bonds). Under pressure, it exhibits H-bond disordering, tunneling, and ﬁnally, H-bond symmetrization at ∼ 18 GPa. This study investigates these 300 K pressure-induced state changes in δ with ab initio calculations. H-bond disordering in δ was modeled using supercell multi-conﬁguration quasiharmonic calculations. We examine: (a) energy barriers for proton jumps, (b) the pressure dependence of phonon frequencies, (c) 300 K compressibility, (d) neutron diﬀraction pattern anomalies, and (e) compare ab initio bond lengths with measured ones. Such thorough and systematic comparisons indicate that: (a) proton “disorder” has a restricted meaning when applied to δ . Nevertheless, H-bonds are disordered between 0 and 8 GPa, and a gradual change in H-bond conﬁguration results in enhanced compressibility. (b) several structural and vibrational anomalies at ∼ 8 GPa are consistent with the disappearance of a particular (HOC-12) H-bond conﬁguration and its change into another one (HOC-11*). (c) between 8–11 GPa, H-bond conﬁguration (HOC-11*) is generally ordered, at least in short- to mid-range scale. (d) between 11.5–18 GPa, H-bond lengths approach a critical value that impedes compression, resulting in decreased compressibility. In this pressure range, especially approaching H-bond symmetrization at ∼ 18 GPa, anharmonicity and tunneling should play an essential role in the proton dynamics. Further simulations accounting for these eﬀects are desirable to clarify the protons’ state in this pressure range.


I. INTRODUCTION
δ-AlOOH (δ) with the CaCl 2 -like or distorted stishovite structure [1] is a critical, high-pressure hydrous phase in subducted slabs in the Earth's lower mantle.The 3D network of strong Al-O-Al ionic bonds reinforces the structure and stabilizes δ in a vast pressure/temperature range.Experiments have shown that δ remains stable up to core-mantle boundary (CMB) pressures, ∼136 GPa [2,3], and temperatures approaching the cold-slab geotherm [4].δ can form solid solutions with other isostructural phases, e.g., MgSiO 4 H 2 phase H [5] and ε-FeOOH [6].The existence of δ substantially extends the stability field of the solid solutions, allowing the ternary system to survive extreme pressures [7] and eventually transport water to the CMB region through slab subduction to participate in the global water cycle [8].
Although δ's crystal structure has a P 2 1 nm space group (No. 30) [1], its H-bond arrangement remained unresolved in powder X-ray diffraction and attracted much interest.Ab initio calculations predicted δ forms asymmetric Hbonds (O-H•••O) at low pressures [9].It also indicated that H-bonds should symmetrize under pressure, forming only ionic (O-H-O) bonds at ∼30 GPa accompanied by an increase in elastic moduli [9].The O-H•••O configura- * rmw2150@columbia.edution consists of a longer hydrogen bond (H-bond), O•••H, and a shorter ionic bond, O-H.Tsuchiya et al. [9] proposed two viable hydrogen off-centered (HOC) unit cell models for δ, HOC-1, and HOC-2.The relative position of two protons in the CaCl 2 -like primitive cell defines these models, with two likely sites for the second proton, H1 or H2, therefore HOC-1 and HOC-2 models.In the H1 site, the second proton relates to the first one by a 2 1 [100] screw operation, while in the H2, it relates by a 2 1 [010].HOC-1 has lower enthalpy in static in ab initio calculations using the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) [10].These calculations predict H-bond symmetrization at 30 GPa.Several experimental studies confirm this asymmetric proton distribution at low pressures [11,12].Still, elastic moduli stiffening signaling H-bond symmetrization is seen at 8 GPa, a much lower pressure than the ab initio static prediction [13][14][15].In addition, proton site occupancies were observed to change under pressure, with H1 sites favored at 0 GPa [12], but H2 sites also being partially occupied at higher pressures [16,17].An order-disordered H-bond transition model was introduced to explain this partial occupancy change with pressure [16,17].The stiffening in elastic moduli at 8 GPa correlates with a redistribution of protons among these sites.Although H-bonds remain asymmetric after this transition, H1 and H2 sites become equally occupied beyond 8 GPa [16].This transition is associated with the extinction of the 0kl neutron diffraction peaks [16,18] of the P 2 1 nm structure.The structure after 8 GPa is referred to as a fully-disordered or the P nnm-disordered model [16][17][18].Further compression would merge the H1 and H2 sites and symmetrize H-bonds at ∼18 GPa [16], which might still be a relatively low pressure compared to the static ab initio (PBE) prediction of 30 GPa [9,19].
H-bond disordering is a general phenomenon in numerous hydrous and nominally anhydrous mineral phases (NAMs) [20][21][22].From a microscopic perspective, H-bond disordering implies the coexistence of multiple H-bond configurations on short-, medium-, or long-range scales.In hydrous phases, this process is usually associated with H-bond symmetrization.As previously shown in H2O-ice [23] and here shown in Secs.III A and III C below, with increasing pressure, there are simultaneous reductions in (a) the energy barrier height between H-sites and (b) in the differences between H-bond configuration energies.Both contribute to the much-increased likelihood of disorder when d(OO) approaches its critical value ∼2.4 Å in both ice-VII [23][24][25] and δ [16].The energy barrier reduction also facilitates proton mobility, especially tunneling, translating into more significant proton redistribution and less well-defined H-sites or d(OH).In addition, proton tunneling has been shown to occur in ice [23][24][25] and is expected to occur in δ [16], although not in the partial H-bond disorder stage (∼5.6 GPa) [26].
The interplay between disordering, tunneling, and Hbond symmetrization, results in a complex multi-stage phenomenon.For example, a series of ab initio studies [23][24][25] have identified the "molecular-ionized-centered" stages throughout the H-bond symmetrization process in the ice-VIII-VII-X system.In the case of δ, although similar stages of the H-bond symmetrization transition have been confirmed by experiments [16] with transition boundaries constrained to ∼9 GPa and ∼18 GPa [16], ab initio calculations have not explored this phenomenon yet.Also, how changes in H-bond arrangements lead to the observed anomalous macroscopic properties such as stiffening in compressibility [14] have not yet been clarified.The pressure evolution of the OH -infrared (IR) and Raman frequencies [13,27] is another issue that needs quantitative interpretation by ab initio phonon calculations.
There is no lack of ab initio studies of δ in recent years.Unfortunately, most studies focus on ordered primitive cells, i.e., HOC-1 and HOC-2 (e.g., see Refs.[19,26,28]), favoring the HOC-1 model because of its lower static enthalpy.Addressing δ with a single configuration disregards the characteristic disorder [16] and its effects entirely.Disorder as the coexistence of multiple H-bond configurations can be modeled by ab initio calculations using supercells.For example, using a 1 × 1 × 2 supercell accounts for the interaction of neighboring H-bond configurations along the c-direction.The adoption of 1 × 1 × 2 supercells [29] allowed the four broad peaks of δ in the OH-stretching frequency regime in the ambient pressure Raman spectrum [30,31] to be successfully modeled.
The success of this simple multi-configuration (mc) model [29] is not accidental.The "fully-disordered" or "P nnm-disorder model" invoked in several experimental studies [16,18] is not entirely accurate either.This is because H-bond disordering is restricted by specific rules, such as "ice rules" in H 2 O-ice [32].Disorder is a complex phenomenon in stoichiometric hydrous phases with high H content because these "ice-disorder"-like rules can vary in each system.In the case of δ, in a regime where proton positions are well-defined, each Al-centered octahedron has a single OH -radical associated with it.As we will show in this paper, both O-H and O•••H bonds remain in the x, y plane.Together, these constraints impose Hbond ordering in two dimensions (a-and b-directions), at least medium-range order (MRO), with disorder limited to a single dimension, the c-direction.Therefore, an ensemble of 1×1×2 supercell configurations might capture the main effects of disorder in δ.The coexistence of H-bond configurations can be treated with the multiconfiguration quasiharmonic approximation (mc-QHA) [33].This approach has been used to address the free energy of ice-VII [33] and has successfully described the order-disorder transition in ice-VIII to -VII in a similar H-bond disordering problem [34].
H-bond symmetrization is convoluted with a spintransition in (Al, Fe)OOH [35] and with Mg-Si site disorder in MgSiO 4 H 2 [36].Because of the abundant and detailed studies of δ's structure at pressures below ∼18 GPa, the H-bond symmetrization pressure, and the unique 1D disorder that allows direct observation of proton distribution projected onto the (001) plane, δ is a convenient case to study this order-disorder phenomenon that also happens in the other related phases, (Al, Fe)OOH and MgSiO 4 H 2 .A better understanding of disordering, tunneling, and H-bond symmetrization in δ would also allow us to more confidently address these more complex cases.
This study seeks to clarify the multi-stage disordertunneling-symmetrization transition process in δ with ab initio calculations.We hope to provide a more detailed atomistic interpretation of macroscopic observations.Although H-bond is itself challenging for ab initio treatments, we would still like to understand the capability and accuracy of the standard DFT functionals for such a problem.For validation against room temperature measurements, we focus on the transition sequence at 300 K. Sec.II first summarizes the methods and calculation details, then tests the accuracy of DFT functionals.Sec.III presents our computational investigation of the multistage process; we address a) energy barriers for protonhopping, b) vibrational properties and the pressure range of validity of the quasiharmonic approximation (QHA), c) order-disorder transition using mc-QHA [33], its implications for the compressibility, and (001)-projected proton distribution d) neutron diffraction anomalies and disorder, and e) bond-length anomalies and tunneling.Sec.IV summarizes our results.

A. DFT Details
Our calculations were performed with the Quantum ESPRESSO (QE) code suite [37] using the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) [10] to density functional theory (DFT).Evolutionary-optimized projector-augmented wave (EPAW) [38] datasets were used to describe valencecore electron interaction for all elements with electronic configurations 3s 2 3p 1 , 2s 2 2p 4 , and 1s 1 being used for Al, O, and H.The plane-wave cutoff energy was set to 80 Ry, and the Monkorst-Pack k-point grid for SCF calculations was set to 4 × 4 × 4 for primitive unit cells and 4 × 4 × 2 for 1 × 1 × 2 supercells.At each pressure, dynamical matrices were obtained using density functional perturbation theory (DFPT) on a 2×2×2 q-point grid and force constants were Fourier interpolated to produce phonon frequencies and vibrational density of states (VDoS) on an 8 × 8 × 8 grid.The climbing-image nudged elastic band method (NEB) [39,40] calculation is performed to determine the energy barrier between different configurations of HOC-δ.
Standard DFT calculations do not easily describe systems with H-bonds.However, δ consists of a 3D interconnected Al-O-Al network formed by strong ionic bonds, which standard DFT functionals, both LDA (local density approximation) [41] and PBE [10], describe well.Although a subtle distortion of the Al-O-Al network is caused by the H-bond [42], the general compressional mechanism of δ is established by the Al-O bond compressibility and should be well described by both these functionals.However, the O-H•••O configurations are more sensitive to the choice of functionals.For example, the H-bond (H•••O) and ionic bond (O-H) lengths are expected to be overestimated and underestimated, respectively, as in ice VII and VIII [34,43], more so with LDA than with PBE.Since this paper primarily addresses the H-bond, we choose to use the PBE-GGA functional.
The validity of the DFT functional can be tested by comparing pressure, i.e., equation of state (EoS), axial lengths, and interatomic distances vs. volume with experimental data.Test calculation results in Figs.S1 and S2 refer to the structures shown in Fig. 1, which will be discussed below.For comparison, LDA results obtained with ultrasoft pseudopotentials [44] are also shown.However, HOC-12 and 11* configurations are unstable with LDA in the volume range of interest and are not shown.There are imaginary phonon frequencies in this volume range.As seen in Fig. S1, static LDA calculations underestimate and PBE overestimate pressure [45], with a difference of ∼5 GPa.Equilibrium axial lengths shown in Fig. S2 S2(f)).Therefore, we prefer the PBE functional to study the H-bond disorder and H-bond symmetrization process of δ-AlOOH.

B. Thermodynamics
Multi-configuration QHA (mc-QHA) calculations were performed with the open-source Python code qha [33].The QHA Helmholtz free energy for each configuration is given by The multi-configuration QHA partition function Z QHA (V, T ) is then obtained by considering multiple inequivalent configurations with free energy F m (V, T ) (Eq. ( 1)) and multiplicity g m [34,46]: The multi-configuration free energy F QHA (V, T ) is calculated from Z QHA (V,T) as The population of the m-th set of configurations, n m (V, T ), is given by [34,46]: Ensemble average of property (A(V, T )), e.g., diffraction intensity I 021 , bond d(OH), are computed as an average of A m weighted by the population n m , i.e., All quantities (e.g., n m , A) can then be converted from a function of (V, T ) to a function of (P, T ) using V (P, T ), i.e., where P (V, T ) = − (∂F QHA /∂V ) T .This method has been successfully applied to obtain the thermodynamic phase boundary between H-bond disordered ice-VII and ice-VIII [34], to study hydrous defects in forsterite [20,46], and to study an order-disorder transition in a low-pressure analog [47] of an ultra-high-pressure planet forming silicate [48].
C. Structure δ-AlOOH's primitive cell contains two formula units (f.u.).Except for results in Sec.III A, this study will model H-bond disorder with 1 × 1 × 2 supercells containing four δ-AlOOH f.u.. Fig. 1 shows four symmetrically distinct 1 × 1 × 2 supercell structures included in this study.These are the same supercells proposed by Tsuchiya et al. [29], except that duplications (i.e., 11*/22* and 12/12*) are merged.The 1 × 1 × 2 supercells can have 16 (2 4 ) configurations in total, which reduce to the four symmetrically distinct configurations shown in Fig. 1 with multiplicities 2, 4, 8, and 2. Any other configuration in larger supercells can be described as a combination of these four sets of configurations in variable amounts.The interaction energy between these configurations is not included in mc-QHA calculations.

D. Pressure notation
At the same volume, PBE/mc-QHA calculations generally overestimate pressure by ∼5 GPa compared to measurements.To avoid confusion, we use the GPa EXP/PBE/QHA unit to refer to pressure from 300 K measurements, static PBE calculations, and 300 K QHA or mc-QHA calculations in subsequent discussions.

III. RESULTS AND DISCUSSION
A. Proton jump energy barrier vs. pressure Before investigating proton disordering in δ, we use NEB to calculate the activation energy required for a single proton jump (or energy barrier) and its pressure dependence.There are two protons in each HOC-1 or HOC-2 unit cell configuration of δ (HOC-11 and HOC-22 supercells in Fig. 1).A single periodic proton jump in the unit cell from the H1 to the H2 site results in a "forward" structural change (from HOC-1 to HOC-2); a jump from H2 to H1 results in a "backward" change.The activation energies for the forward and backward periodic jumps in a single unit cell are shown in Fig. 2(a) as solid and dashed blue curves, respectively.At 0 GPa PBE , the activation energy is ∼0.075 eV (∼870 K).Under compression, the energy barrier decreases significantly to roughly 0.025 eV (∼300 K) at ∼11 GPa PBE and almost vanishes at ∼20 GPa PBE .The barrier height reduction justifies the evolution of proton distribution in a classical Langevin molecular dynamic (MD) analysis that shows a change in OH radial pair-correlation distribution (Fig. 4 in Ref. [49]) from discrete at 0 GPa PBE to interconnected at ∼10 GPa PBE [49].
A path integral MD analysis of proton motion shows that protons in H-bonds would be able to tunnel through the barrier wall roughly when the barrier height between neighboring proton sites become comparable to k B T , the thermal energy [23].Our calculations show that at 0 GPa PBE , the energy barrier for a forward jump is larger than ∼0.026 eV (∼300 K).But at higher pressure, the barrier decreases to less than k B T at 300 K, showing that δ could enter a strongly anharmonic and tunneling regime before or around 11 GPa PBE .From a statistical perspective, p jump /p base = exp(−∆E/k B T ) where ∆E = E jump − E base ).At 300 K and 0 GPa PBE , ∆E ∼ 0.08 eV and p jump /p base ∼ 1/25; at 11 GPa PBE , p jump /p base ∼ 1/2.7 only.At higher pressures and temperature (P, T ) conditions of subducted slabs [4], it is entirely possible for the structure to disorder (as described below).To simulate a solitary proton jump in a crystal more accurately and understand size effects, we model a single proton jump in various supercells.They are 1 × 1 × 2, 1 × 1 × 4, and 1 × 2 × 2 multiples of a unit cell.Fig. 2(b) represents the "base" and "jumped" configurations.For one jump on a 1×1×2 supercell, the initial and final structures correspond to HOC-11 and HOC-12 in Fig. 1.The energy barrier increases for a forward jump and decreases for a backward jump.At ∼17 GPa PBE , the backward energy barrier decays to zero.Other 1 × 1 × 2 supercell final configurations are possible but correspond to multiple jumps.For example, the transition from HOC-11 to HOC-11* (-11, -11*, etc., subsequently) structures in Fig. 1 can be achieved by jumping twice, -11→-12→-11*.
Further increase in supercell size along the c direction does not result in substantial differences in energy barrier (e.g., in 1 × 1 × 2 or 1 × 1 × 3, etc).
Proton jumps on other types of supercells offer surprising results.We are unable to find a stable static structure with a single jump in (001), the a, b-plane, e.g., in 1 × 2 × 1 or 2 × 1 × 1 supercells.The 1 × 2 × 2 supercell can produce a single jump final structure at ∼0 GPa PBE , but only in a narrow pressure range (see Fig. 2(a)).Even in this narrow range, the energy barrier for backward jump decreases substantially compared to the other 1 × 1 × N jumped supercell configurations.This result means an isolated proton jump in (001) is unlikely to result in a long-lived state.This result confirms that a single proton jump in an ordered structure is not likely, and proton jumps must be coordinated in the a, b-plane.Even in H-disordered δ, periodicities along the a and b axes are likely to be maintained in some length scale.For this reason, we must consider a more orderly type of disorder.Dynamically, the change in H-bond configuration (disorder) in the a, bplane might occur in sequential jumps as in-plane soliton propagation, which is a consequence of the "one proton per octahedron" limit imposed by the chemistry.
Because the size effect is sufficiently mitigated with supercells with two units along [001], and periodicity should remain mostly intact in the (001) plane, at least within a finite length scale, our results suggest that 1 × 1 × 2 supercell calculations offer ample base to address structural disorder in δ.Moreover, since Ref. [29] has confirmed that the Raman spectrum in the OH-stretching region is generally reproduced with these units, the ensemble of all possible 1 × 1 × 2 supercells might suffice to describe all possible short-range O-H•••O configurations.The six structures given in Ref. [29] can be reduced to four independent ones shown in Fig. 1, HOC-11, HOC-11*/22*, HOC-12, and HOC-22 (11*/22* and 12/12* reported in Ref. [29] are symmetrically related).The remaining sections of this study use these supercells to investigate the compressional behavior of H-bonds in δ.

B. Pressure dependence of VDoS, proton stability, and tunneling signature
After confirming that these supercell models describe well the predominant short-range proton configurations, it is essential to clarify the pressure dependence of their related phonon frequencies and modes.The high energy barriers discussed in the previous section suggest minor anharmonic effects on proton motion at low P, T s.However, these modes should become anharmonic with increasing P, T s.Therefore, the highest-pressure results reported in this section have more phenomenological and qualitative significance than predictive power.Nevertheless, they offer a detailed view of trends.
Fig. 3 shows the evolution of various phonon properties under pressure.Because phonon frequency behavior is a function of volume or bond-length (see, e.g., [43]), and our angle between H-displacement vector and OH bond direction vs. volume.Symbols are from PBE results, and lines are spline interpolations to help identify mode-continuity with volume.The main bottom axis indicates the unit cell volume (2 f.u.); Extra axes at the top show the experimental and PBE static pressure at the investigated volumes.The additional axis at the bottom, d(OO) from static PBE calculations, is given to facilitate discussions.Scattered "+" and "×" symbols in panels (a-h) denote the frequency of IR and Raman peaks extracted from Refs.[27] and [13].
system of interest, O-H•••O, is relatively localized and could be studied vs. d(OO) bond lengths [23], it is more informative to display frequencies vs. volume henceforth.As will be shown in Sec.III E, d(OO) is nearly a linear function of volume and is well-predicted.Therefore, we add a scale showing the corresponding PBE d(OO) bondlengths.We also supply different scales at the top, indicating PBE and experimental pressures at corresponding volumes.Phonon frequency continuity between adjacent volume points is established using eigenvectors similarity described in Ref. [50].
Figs. 3(a-d) show the pressure evolution of the vibra-tional density of states (VDoS) obtained by first interpolating phonon frequency vs. volume, then calculating the VDoS at each volume on a dense volume grid for each of the four structures in Fig. 1.Across the spectrum, the phonon bands cover three broad frequency domains: the 500-800 cm −1 bands consist of Al and O (heavier-ions) related modes, the 1000-1500 cm −1 bands are OH bending modes, and the 1500-2500 cm −1 ones are OH stretching modes.The pressure dependence of the 500-800 cm −1 modes is similar to those of stable anhydrous phases and overlaps with Raman mode frequencies [13].The high-frequency modes, both the 1000-1500 cm −1 and the 1500-2500 cm −1 bands, are more strongly pressuredependent.The stretching and bending modes' negative and positive pressure coefficients are similar to those reported in previous theoretical studies of δ [29,49] and other H-abundant phases [51,52].The current static calculation slightly underestimates d(OO) and the stretchingmode frequencies compared to other works [29,49].All phonon frequencies exhibit a near-linear dependence on volume until the OH-stretching and -bending mode frequencies become similar and these modes start interacting.
Comparing the VDoSs of the four supercell models reveals a meaningful relationship between the OH stretching band frequencies and proton arrangements along the interstitial channels.There are two stretching mode bands, one beginning above and one below 2000 cm −1 at 0 GPa PBE .The former exists in the HOC-11, -22, and -12 configurations and the latter exists in HOC-11* and -12 (see Fig. 3(a-d)).There are two types of proton arrangements: in HOC-11 and -22, protons are aligned along the interstitial channels parallel to the c axis, i.e., their positions are related by a translation c.In HOC-11*, protons are half-aligned, i.e., their positions are related by a 2c translation.This proton arrangement is more "disordered".The HOC-12 configuration has both aligned and half-aligned types of protons.Aligned protons have higher OH-stretching frequencies than half-aligned ones.This interpretation sheds light on the origin of results obtained with a 300 K quantum thermal bath (QTB) MD simulation performed on 2 × 2 × 4 supercells [49].Two series of OH-stretching modes in the 0-10 GPa pressure range were reported indicating the coexistence of aligned and half-aligned arrangements.The disappearance of the higher-frequency band at 10 GPa can be attributed to the disappearance of the aligned proton configurations, indicative of growing disorder in the structure.Mostly HOC-11*-like arrangements are expected to be present beyond 10 GPa in those simulations.This configuration change will be addressed with mc-QHA calculations in Sec.III C.
Figs. 3(e-h) show the evolution of OH zone center (Γpoint) mode frequencies under compression.Red, green, and blue denote the OH-stretching, OH-bending in the a, b-plane (in-plane bending), and OH-bending along the c-direction (out-of-plane bending).This classification of proton modes is based on local mutually perpendicular axes defined by (a) the in-plane ionic O-H bond direction, (b) the perpendicular direction in (001), the a, b-plane, and (c) a third out-of-plane direction perpendicular to both, i.e., along [001].The intensity of color represents the magnitude of eigenmode projections onto these OH displacement directions.Translucent curves mean the proton displacement in those modes correlates strongly with Al or O ionic motions.
The four proton configurations are stable in different pressure ranges in static ab initio calculations.When a proton configuration becomes unstable, it transforms into another one.Imaginary frequencies in Figs.3(e-h) identify proton vibrational instabilities at the Γ point.HOC-12, -11*, and -22 develop signs of vibrational instability at ∼12, 22, and 26 GPa PBE , or ∼8, 12, and 20 GPa EXP , while HOC-11 survives beyond 26 GPa PBE .Figs. 3(e-h) show that unstable modes have a slightly reddish color, indicating the onset of OH tunneling behavior (stretching mode instability), similar to earlier reports [49].
As for the high-frequency modes, we observe a "mode mixing" and "frequency crossing" between stretching and in-plane OH-bending modes preceding the vibrational instabilities in all four supercells shown in Figs.3(e-h).As stretching mode frequencies approach the in-plane bending frequencies, they repel each other, these modes mix, and the stretching and bending branches start to switch order.This process is continuous, and the eigenmodes neither stretch nor bend throughout a few GPa while these frequencies overlap.Similar mixing behavior was previously reported in these modes' harmonic analysis [29,49].An OH-stretching bandwidth increase (see Figs. 3  (a-d)) accompanies this "mode mixing".
This mode "mixing" process is captured in the volume dependence of the angle between the O-H displacement vector and the O-H bond direction shown in Figs.3(i-l).Fig. S3 shows in more detail the pressure evolution of these eigenmodes in HOC-11*.Over the entire process, the displacement directions of both modes undergo a continuous 90°rotation.We do not observe any strong distortion in the O-octahedra throughout this process.This switch in the mode frequency order results from the contraction of interstitial site volume or H-bond length under pressure.A similar phenomenon was reported in MD simulations of disordered ice-VII at d(OO) around 2.4 Å [23], and it could be a general phenomenon in O-H•••O bearing systems.In ice-VII [23], this mode crossing behavior was described as a pre-tunneling effect, suggesting that δ with specific H-bond configurations will enter a similar regime at corresponding pressure (e.g., 11 GPa PBE or 8 GPa EXP for HOC-12, 16 GPa PBE or 10.5 GPa EXP for HOC-11*).We will address this issue further in the upcoming sections under the context of proton configuration population evolution from the mc-QHA analysis.
With a better understanding of the pressure dependence of eigenmodes and their frequencies, we can reinterpret the pressure dependence of IR frequencies [27] and tentatively connect them with disordering, tunneling, and symmetrization in subsequent sections.When the measured IR OH-bending mode frequencies [27] shown in Fig. 3 are more clearly displayed in Fig. S4, various regions with different frequency vs. volume dependencies are better identified.Slope changes indicate changes in the proton arrangement or vibrational state in δ.For in-plane OH-bending frequencies, the first slope change happens at ∼7.5-8 GPa EXP .There is an abrupt compressibility increase [14], and the 021 neutron diffraction peak disappears [16,18] at this same pressure.It also happens that the HOC-12 configuration becomes vibrationally unstable at this pressure (see Figs. 3(c,g,k)).The disappearance of the 2800 cm −1 branch at 10 GPa in MD simulations [49] should also correspond to this change because the HOC-12 configuration essentially disappears after this point.The slope change in the out-of-plane OH-bending frequency happens later, at ∼11 GPa EXP (see Fig. S4(a)).A kink is observed immediately after at 11.5 GPa EXP , marking the beginning of the 11.5-18 GPa EXP range in which d(OH) distances behave anomalously [16] (see Fig. S4(b)).We correlate this change with the anharmonicity in the HOC-11* type structure which starts at ∼18 GPa PBE (∼11.5 GPa EXP ).Finally, at 18 GPa EXP , we observe slope changes in both in-plane and out-of-plane bending frequencies (see Fig. S4).This pressure corresponds to the full H-bond symmetrization pressure [16].These connections will be explored in greater detail in subsequent sections and give fresh insights into the nature of this complex multi-stage phase change in δ complicated by proton tunneling.
C. Population analysis, thermodynamic properties from mc-QHA, 300 K EoS, implications for compressibility, and proton distribution The small (< 0.1 eV/f.u.) energy differences between the four supercell H-bond configurations shown in Fig. 1 suggest the disordered system is indeed possible [29].We now investigate the thermodynamic properties of HOC-δ using the mc-QHA method [33] and the volume dependence of VDoSs of the four supercells.As previously discussed in Sec.III A, the four structures shown in Fig. 1 represent well possible short-range proton configurations.Their coexistence in different proportions describes structural disorder in our model.
Because QHA calculations require stable phonons, four separate mc-QHA calculations were performed in different pressure ranges (0-12, 12-20, 20-26, and 26-30 GPa QHA ).Different proton configurations are stable in these pressure ranges.Since phonons are calculated assuming harmonic potentials, anharmonic effects and proton tunneling are not accounted for in this treatment.The current mc-QHA approach can shed light on disordering effects and H-bond symmetrization but is not predictive.
In this section, thermodynamic properties are discussed at 300 K and QHA pressure.However, our 300 K QHA calculation systematically overestimates pressure by ∼4-6 GPa compared to experiments at the same volume due to the use of the PBE exchange-correlation functional and other factors, e.g., anharmonicity (see Fig. 4(c) and Fig. S1).Therefore, a shift in pressure is needed when comparing our proposed boundary pressures with experimental pressures (e.g., −6 GPa shift applied to Fig. 4(c)).This pressure shift is also required when comparing measured and calculated d(OO) bond lengths.
The vibrational instability in the HOC-12 structure at ∼12 GPa QHA leads to the disappearance of this atomic arrangement beyond this pressure.In the 12-20 GPa QHA range, HOC-11* becomes the most abundant configuration, and the HOC-11's and HOC-22's populations remain low.After HOC-11* becomes unstable at ∼20 GPa QHA , HOC-22 and -11 configurations are still stable in the 18-26 GPa QHA range.In the 26-30 GPa QHA range, near H-bond symmetrization, only HOC-11 seems to survive, but HOC-22 and HOC-11 are the same when the H-bond is symmetric.
A closer look at the free energy of these structures (Fig. S5) sheds light on the driving force behind this struc-tural evolution.We showed in Sec.III B that H-bond arrangements affect stretching mode frequencies.Fig. S5 shows how these differences in vibrational properties contribute to the free energy, in particular, the zero-point motion energy (E zp ).Because H-bond arrangements halfaligned along the [001] have lower stretching mode frequencies, half-aligned proton configurations have lower E zp .Although HOC-11 has the lowest static energy/enthalpy, the vibrational energy contribution becomes more important as differences in static energy decrease under compression.As shown in Fig. S5, from ∼5 GPa PBE onwards (or ∼1 GPa EXP ), HOC-11* and -11 become the most and least stable configurations, respectively, owing to E zp .With mixed proton configurations, HOC-12 is not at the extremes of this energy spectrum.Still, its greater multiplicity makes it the most populous configuration in the 0-12 GPa EXP range until it becomes unstable at 12 GPa EXP .After 12 GPa EXP , the HOC-11* configuration is the most abundant due to its energy and multiplicity.Essentially, the decreasing static freeenergy difference in the disordered system allows the half-aligned proton configurations to be more populous at higher pressures after the HOC-12 configuration becomes unstable.In summary, disregarding anharmonicity or dynamic disorder, the current mc-QHA calculation suggests the transitions in HOC-δ occur in the following order: in 0-8 GPa EXP (0-12 GPa QHA ) δ undergoes a continuous change in H-bond configuration with HOC-11 transforming into -11*; at 8 GPa EXP (12 GPa QHA ), HOC-12 transforms to -11*; at 14 GPa EXP (20 GPa QHA ), HOC-11* transforms to -22; at 20 GPa EXP (26 GPa QHA ), HOC-22 transforms to -11.Because of multiplicity and zero-point energy, the 0 GPa EXP measured Raman spectrum [31] resembles the calculated HOC-12 spectrum more closely than the -11 [29], the most stable configuration at this pressure.This underlines the necessity of including configuration multiplicity and vibrational contributions in ab initio calculations of H-bond disordered systems.
The retirement of the aligned protons in the HOC-12 configuration at 8 GPa EXP , i.e., the instability in HOC-12 and the dominance of the -11* configuration afterward, is consistent with the disappearance of high-frequency stretching bands at ∼10 GPa in 300 K QTB MD simulations [49].At least within a limited pressure range, i.e., up to 12 GPa EXP (18 GPa QHA ), the predominance of the half-aligned proton configuration in HOC-11* is confirmed by both mc-QHA and QTB MD simulations at 300 K [49].After 12 GPa EXP (18 GPa QHA ), mc-QHA calculations and QTB simulations at 300 K [49] no longer offer consistent results, most likely because mc-QHA calculations cannot properly capture anharmonic effects.In the QTB MD simulations [49], no discontinuity is observed in the phonon stretching band, which is expected if the half-aligned proton configuration in HOC-11* vanishes and the proton-aligned configurations in HOC-11 and -22 become dominant as predicted by mc-QHA.Despite the softening of branches approaching imaginary frequencies (Fig. 3(f)), we can still optimize the HOC-11* structure or a distorted version in static calculations.This signals a small energy barrier (Fig. 2(a)) that, in reality, might fail to trap protons locally due to tunneling or anharmonic fluctuations.
This population evolution vs. pressure significantly affects the 300 K EoS and compressibility of the multiconfiguration system, as seen in Fig. 4(b).This figure also shows the 300 K compression curves for each of the four configurations for comparison.Their volumes differ by less than 0.4 Å3 at any pressure and are in the following order: V 11 > V 22 > V 12 > V 11 * .Intuitively, the half-aligned H-bond configurations keep protons farther apart, lower their stretching mode frequencies, and reduce the volume.
Piecewise, the overall compression curve predicted by mc-QHA resembles more closely that of the predominant configuration in different pressure ranges (see Fig. 4(c)).In the 0-12 GPa QHA (0-9 GPa EXP ) range, the decrease in population of aligned configurations (HOC-11 and -22) and the increase of half-aligned ones (HOC-12 and -11*) results in an increase in compressibility w.r.t. that of any of the configurations alone (see Fig. 4(b)).The same happens with the extinction of HOC-12 at 12 GPa QHA .This result correlates well with the measured high compressibility regime (0-9 GPa EXP ) before the anomalous stiffening [14].
In a limited pressure range of 12-20 GPa QHA or 9-14 GPa EXP , HOC-11* alone determines the shape of the compression curve.In this pressure range, the compressibility decreases compared to that in the 0-12 GPa QHA range (see Fig. 4(b)).This behavior correlates with the observed behavior in the 9-12 GPa EXP pressure range [14].
In the 20-26 GPa QHA stage (or ∼12 GPa EXP onwards), mc-QHA predicts HOC-22 and -11 configurations determine the overall compressibility.These configurations have larger volumes.The ∼0.4 Å3 discontinuous volume increases at 20 and 26 GPa QHA (see Fig. 4(c)) should somehow be smoothed in the QHA compression curve leading to a further decrease of compressibility, as it seems to occur at 14 GPa EXP in experiments [14].More advanced quantum simulations are desirable to shed light on proton states and configurations in this pressure range.This will be discussed in greater detail when the H-bond's compressional behavior is addressed in Sec.III D.
Although PBE calculations systematically overestimate pressure by ∼4-6 GPa or volume by ∼2 Å3 , our overall segmented mc-compression curve shape in Fig. 4(c) seems consistent with the experimental compression curve [14,16,18], particularly in the 0-14 GPa QHA range.Because protons should be tunneling at and beyond ∼18 GPa QHA (12 GPa EXP ), our results are not predictive beyond this pressure.However, the phase consisting of HOC-11*, HOC-11, and HOC-22 components above ∼15 GPa QHA displays an overall compression behavior in good agreement with experiments.

Implications for the (001)-projected proton distribution
As indicated in Sec.III A, unconstrained disorder is unlikely in individual a, b-planes, at least in short-to medium-range scale.Disorder is thus restricted to the cdirection.Previously measured proton distribution maps vs. pressure (Fig. 3 in [16]) correspond to projections of the proton distribution on (001).Here we examine the correlation between the (001)-projection of the proton arrangements to the state of disorder along the c-direction.
Consider single crystals of each of the four HOC structures in Fig. 1.The HOC-11's and -22's (001)-projections will show a single proton distribution peak.HOC-11*'s will show two symmetric proton distribution peaks; HOC-12's will have two asymmetric peaks.Next, consider entirely random stackings of equivalent configurations for each of these HOC units along [001].The (001)-projected proton distributions will show two symmetric proton distribution peaks for any of these randomly stacked configurations.These are, of course, the two extreme cases.The coexistence of multiple configurations would introduce an additional degree of disorder.Therefore, asymmetric peaks in the (001)-projected proton distribution map require, at least to some extent, preservation of some order along [001].
"Partial disorder" as measured at 6.37 GPa EXP (see Fig. 3(a) in [16]) best describes the proton distribution at low pressures (i.e., < 8 GPa EXP ).The proton sites in the 6.37 GPa distribution map are discrete, showing H-bonds are well-defined.Some ordering in the predominant configuration, i.e., HOC-12, could also explain the asymmetric peak distribution.The continuous population change from HOC-11 to -11* and the increasing disordering trend could reduce peak asymmetry [16] under compression.A recent low-field 1 H-NMR measurement at 5.6 GPa EXP [26] confirms that H-bond tunneling does not occur in this pressure range, leaving H-bond disorder the predominant effect at low pressure.
Around 8.4 GPa EXP and below 9.5 GPa EXP (Fig. 3(b) in [16]), the proton distribution peaks show broad and asymmetric distribution [16].This behavior signifies that protons are no longer localized at one site but can move through/across the energy barrier, which should be low (see Sec. III A) and have a more gentle curvature.Based on the pressure of this anomaly and the disappearance of OH-stretching band in 300 K QTB MD [49] at corresponding pressure, we correlate this behavior to HOC-12's instability (depicted in Fig. 3(g)) and proton's tunnelinglike behavior in a narrow pressure range preceding this instability (Sec.III A).Therefore, the 8.4 GPa EXP (001)projected distribution map seems to have captured an intermediate state which serves as pathway for the HOC-12 to HOC-11* (aligned to half-aligned) configuration change.
A bimodal distribution with the same peak height at slightly higher pressure is observed at 9.5 GPa EXP [16].This behavior is consistent with the predominance of the HOC-11* configuration starting at ∼11 GPa QHA as indicated by the mc-supercell calculation.It is also compatible with the absence of a stretching band in the QTB MD simulation [49].
The current study cannot fully address δ's state after HOC-11* shows signs of tunneling and instability (starting at ∼18-20 GPa QHA or ∼12-14 GPa EXP ).The increase in OH bond-length (next section) and the lower compressibility after -11*'s instability seem to suggest that a transformation back into HOC-11 and -22, the proton-aligned configurations, is possible.However, the existence of these configurations is inconsistent with the absence of OH-stretching band frequencies in 300 K QTB MD [49], which suggests an everlasting dominance of the half-aligned configuration up to the H-bond symmetrization pressure.This behavior could be caused by strong anharmonic effects and tunneling not included in our study.Nevertheless, the dominance of the half-aligned configuration (HOC-11*), disordered mixtures of aligned ones (HOC-11 and -22), or a dynamically disordered δ could all contribute to the (001)-projected proton distribution showing two symmetric peaks before full H-bond symmetrization.
Suppose our theory represents approximately the behavior of proton configurations under pressure.Our results agree that δ is partially disordered at low pressures [16,18].However, our results suggest that the previously proposed "fully-disordered" stage [16,18] might, in fact, be a predominantly half-aligned HOC-11*-type configuration starting at around 8.5 GPa EXP .

D. Evolution of neutron diffraction intensity and the 9 GPa transition
Previous studies [16,18] constrain the pressure of the partially disordered to fully disordered transition to ∼9 GPa EXP by tracking the pressure dependence of the 021 neutron diffraction peak intensity.Because this peak exists in the "partially-disordered" model and not in the "fully-disordered" model, the evolution of the 021 peak intensity should track the state of disorder.However, the configuration population analysis in Sec.III C suggested the 9 GPa EXP configuration change could be better explained by the disappearance of HOC-12-like proton configurations and its transformation into the HOC-11*-like configurations ("HOC-like configurations" means domains, large or small, of these ordered configurations).For consistency, we need to demonstrate a change from partially disordered HOC-configurations to half-aligned configurations could also result in the same peak extinction behavior.[16].The relative 021 peak intensities for HOC-11* and -22 are less than 0.002 and, therefore, not plotted.The "overall" peak intensity is the average of peak intensities from different configurations using the mc-QHA populations at 300 K. Diffraction patterns, peak positions, and peak intensities were calculated using CrystalDiffract® [53] and the optimized crystal structures.
A structure factor analysis can determine the extinction condition for the 021 diffraction peak (F 021 , or F 022 for the 1×1×2 supercell [54]) of the four HOC configurations.Results show that F 021 = 0 for HOC-11 (primitive cell space group P 2 1 nm, No. 31) and HOC-12 (supercell space group P m, No. 31).F 021 = 0 for HOC-22 (primitive-cell space group P n21m, No. 31) and HOC-11* (supercell space group P 2 1 2 1 2 1 , No. 19).In the latter structures, the corresponding peak is absent.Table S1 gives detailed information on proton site symmetries in these structures.With the help of CrystalDiffract® [53], we create the diffraction profiles for ab initio optimized HOC structures and get d-spacings (d hkl ) of the diffraction peaks.Fig. 5(a-d) shows the normalized diffraction patterns of the four supercells at 0 GPa PBE .Most features in our diffraction spectra agree with those in neutron diffraction patterns [16].An observed systematic trend is that low d-spacing peaks have low diffraction intensities [16].This trend is not present in our calculated diffraction patterns, but it should be irrelevant to our analysis focused on the pressure evolution of one particular peak intensity.The calculated d 021 spacing is ∼1.7 Å, consistent with that reported [16,17].In Fig. 5(a-d), d 021 's are indicated with arrows.d 021 's differ slightly for different HOC structures at the same PBE pressure because lattice parameters differ slightly.CrystalDiffract® [53]'s prediction shows that the 021 peaks exist only in HOC-11 and -12 but not in -11* and -22, consistent with our structure factor analysis.Therefore, transformations from HOC-11 or -12 configurations into HOC-11* will extinguish the 021 diffraction peak.
To further illustrate our analysis and to compare it with the reported compressional behavior of the 021 peak [16], we plot the peak intensities ratio, I 021 /I 110 , vs. volume in Fig. 5(e).Our results show I 110 's are independent from H-bond configurations, therefore it is reasonable to use it as the standard for normalization.The I 021 /I 110 ratio is smaller than 0.002 for HOC-11* and -22 over the entire pressure range, thus not shown in Fig. 5(e).The I 021 /I 110 ratios for both HOC-11 and -12 decay linearly under compression, but neither is absent at ∼8 GPa EXP .For this reason, we confirm that a configuration change, continuous or discrete, from HOC-11 or -12 is likely at 9 GPa EXP .As shown in Sec.III C, it is the HOC-12 configuration that disappears at this pressure.The HOC-11 population decreases in favor of the HOC-11*'s up to 9 GPa EXP , but it survives beyond this pressure and is the last one to disappear in our calculations.
Considering the contribution of each configuration based on their 300 K populations (see Fig. 4(a)), a curve reflecting our estimated population-averaged relative intensity ratio I 021 /I 110 is also plotted in Fig. 5(e).Compared with Ref. [16], our estimation of I 021 /I 110 in the 0-9 GPa EXP range is slightly smaller than the experimental measurement, but the pressure dependence of this average I 021 /I 110 shows the correct decaying trend up to 9 GPa EXP .The disappearance of HOC-12 at ∼9 GPa EXP and the dominance of HOC-11* afterwards makes I 021 /I 110 decay to below 0.002 or become effectively absent above 9 GPa EXP , in agreement with Ref. [16]'s observation.The decay and extinction of I 021 /I 110 predicted by the mc-QHA calculations is clearly associated with changes in configuration populations and the disappearance of the HOC-12 configuration at ∼9 GPa EXP (∼13 GPa QHA ).Therefore, it supports the partially disordered to predominantly HOC-11*-ordered proton configuration change we described in the Sec.III C.However, the present study is unable to describe the exact configuration change behavior near the vibrational stability limit, therefore only a sharp transition is shown in Fig. 5(e) instead of more gradual I 021 /I 110 ratio decay as measurements [16,17].As described in Sec.III C, the transition from HOC-12 to -11* should undergo a brief regime of H-bond tunneling, where protons have a large position spread [16], weakening the I 021 diffraction intensity as a result.At this point, tunneling might not have developed to its full-fledged behavior before HOC-12 transforms into -11*.Tunneling here should be viewed as a pathway for change in H-bond arrangement only.Further experimental and theoretical investigations are needed to detail this proton configuration change phenomenon at 9 GPa EXP .In addition to other evidence presented earlier, the good agreement between the predicted and measured development of the I 021 diffraction peak intensity provides further evidence that the observed ∼9 GPa EXP change could correspond primarily to a change from predominant HOC-12 to -11* configuration.
After a brief regime of HOC-11* stability, mc-supercell calculations predict HOC-11* transforms into HOC-22 and -11 before full H-bond symmetrization.If so, the reappearance of the 021 peak might be expected alongside -11's reappearance at ∼26 GPa QHA (∼18 GPa EXP ).This effect has not been reported.Suppose tunneling or anharmonic fluctuations are happening between HOC-22 and -11.In that case, H-bond may appear symmetric, as implied by merging the two symmetric H-distribution peaks in the 18.1 GPa EXP (001)-projection map [16].This behavior is actually suggested by the 300 K QTB MD showing a single high-frequency H-bond stretching bandwidth [49].In this case, the 021 neutron peak is not expected to reappear.

E. Interatomic distances d (OH) and d (OO) and tunneling
This section investigates the evolution of interatomic distances d(OO) and d(OH) under pressure.Fig. 6 shows the volume dependence of d(OO) and d(OH) bond-lengths.We first notice that d(OO) depends almost linearly on volume in the 0-20 GPa EXP pressure range.Whether measured [16] or computed, this bond-length is undisrupted either by protons' anomalous behavior [16] or by the differences between H-bond arrangements in ab initio Between 9.5-11.5GPa EXP (∼15-18 GPa QHA ), , protons are not tunneling because the 9.5 GPa EXP measured (001)-projected proton distribution shows a welldefined double-peak shape [16] (b,f,j)), structural instability in ab initio calculations, and discontinuities in measured bond-length trends [16].Therefore, the ∼11 GPa EXP discontinuities in measured bond-lengths [16] and the ab initio instability in HOC-11*, possibly accompanied by a pre-instability tunneling behavior, are likely correlated.If so, our static ab initio calculation predicts well the d(OO) and both d(OH)s limiting values.This d(O-H) limit is reflected on the measured OH-bending mode frequency [27] vs. volume.This limit separates two regimes distinguished by this mode frequency's pressure/volume dependence (see Fig. S4(a)).In ab initio phonon calculations, OHstretching and in-plane bending bands start mixing at this point (see Sec. III B, Figs.3(e-h)).

F. Implication
Our work shows that the H-bond disordering process that precedes H-bond symmetrization can be described as a sequence of changes in predominant H-bond configurations.Subtle changes in the atomic structure under pressure, as found here, contribute to anelastic effects.We identify such an effect in the pressure range of ∼10-20 GPa QHA , which corresponds approximately to the pressure range of ∼6-15 GPa EXP .In the latter pressure range, Brillouin scattering experiments identify a rapid increase in acoustic velocities and anomalies in the compression curve.Here we suggest that this rapid increase in acoustic velocities is related to anelastic effects, i.e., structural accommodation under pressure, before H-bond symmetrization.Anelastic effects are known to cause seismic wave attenuation.As previously suggested [13], fast velocity anomalies in the transition zone (15-25 GPa EXP ) may signal the presence of δ in this region.Here we suggest that if seismic attenuation is observed concurrently with fast seismic velocities (low-frequency acoustic waves) in the transition zone (410-660 km depth), it will strengthen the evidence for the presence of δ in this region.High crust/slab temperatures in Earth's interior might further enhance H-bond disordering, and stiffening in compressibility will likely occur at even lower pressures.The magnitude of these effects in the mantle depends on the relative abundance of δ, which correlates with the degree of slab hydration and is quite uncertain.
These predicted anelastic effects are expected in other structurally similar water carriers in the mantle, e.g., MgSiO 4 H 2 phase H [5], -FeOOH [6], and Al,H-bearing stishovite.They might also be expected in other dissimilar phases since the phenomena addressed in this paper are typical of H-bonds.The importance of these effects to seismic velocities will be proportional to the abundance of the hydrous phase in the hydrated slab, which is uncertain [55].Therefore, a precise characterization of δ-AlOOH's acoustic velocities and other co-existing hydrous phases will also help clarify the degree of slab hydration in the mantle up to CMB depths.
At 300 K, our calculations are not predictive beyond 18 GPa QHA (∼12 GPa EXP ), where anharmonic effects become strong.Therefore, it is useless to speculate further on their geophysical consequences.Predictive calculations beyond this pressure will require genuinely quantum MD simulations to address H-bonds at high temperatures, encompassing the superionic and the dehydration regimes [2,56] that occur near ambient mantle temperatures.

IV. CONCLUSION
In summary, using multi-configuration quasiharmonic (mc-QHA) ab initio calculations, we studied a sequence of H-bond configuration changes in H-off-center δ-AlOOH (HOC-δ) under pressure at 300 K, up to full H-bond symmetrization.
The energy barrier for proton jump calculations confirms the notion that H-bond disordering is limited by the presence of one OH -per aluminum octahedron.This constraint imposes short to medium range H-bond ordering in the a, b-plane and limits disorder to the c-direction.Short-range interaction between a, b-planes with different H-bond arrangements allows us to use a 1 × 1 × 2 supercells to model the disordered δ phase.Therefore, we are able to use four symmetrically inequivalent supercell configurations (16 configurations in total), HOC-11, -12, -11*, and -22, to investigate the change in configuration population vs. pressure with mc-QHA.HOC-11 and -22 contain protons fully aligned in the [001] interstitial channels, HOC-11* has only half-aligned protons, while HOC-12 has both proton arrangements (see Fig. 1).This model consistently reproduces a significant number of experimental observations and allows us to correlate some measured properties with underlying atomistic processes.We do this by inspecting calculated and measured properties vs. volume rather than pressure.
Below 8 GPa EXP , there is considerable variation in configuration populations (see Fig. 4(a)), but HOC-12, with two kinds of proton alignments, is the most abundant one.mc-QHA calculations reproduce well both d(O•••H) and d(O-H) bond lengths (see Fig. 6) and δ's more compressive behavior in this pressure range [14] (see Figs. 4(b,c) and S4(d)).They result from the continuous change in H-bond configuration populations under pressure.Also, the coexistence of multiple configurations (mc) confirms δ's partially disordered nature, which should result in the asymmetric proton distribution projected in (001) [16].mc coexistence also results in the four broad peaks observed in 0 GPa Raman spectra [29,31] and the two OHstretching bands observed in 300 K QTB-MD simulations [49], one from aligned and the other from half-aligned proton configurations.
At ∼8 GPa EXP , the 021 neutron diffraction peak disappears [16] (see Fig. 5(e)), the proton distribution peaks become ill-defined [16], and the pressure dependence of the out-of-plane OH-bending mode frequency changes [27] (see Fig. S4(a)).Also, the higher OH-stretching frequency band disappears in the 300 K QTB-MD simulations around ∼10 GPa [49].We correlate these observations with the vibrational instability and disappearance of the predominant HOC-12 configuration (see Figs. 3(c,g,k), 4(a), 5(e), and S4(a)) and its replacement by -11* with only half-aligned protons.The obscure proton distribution peaks might result from a brief proton tunneling behavior preceding this configuration change.
Between ∼8 and ∼11.5 GPa EXP , proton distribution maps show two symmetric discrete peaks that have been attributed to a fully disordered proton configuration [16], while 300 K QTB-MD simulations display a single OHstretching frequency band [49].Overall, δ becomes less compressible in this brief pressure range [13,14] 6).We correlate these observations with the pre-dominance of HOC-11*, small abundances of -11 and -22 in this pressure range, and the eminent instability of HOC-11* at ∼11.5 GPa EXP .HOC-11*, the densest configuration (see Fig. 4(b)), produces two discrete symmetric peaks (see Fig. 1), is less compressible, and displays smaller compressibility, and a single and lower OHstretching frequency band (compare Figs. 3(b) and 3(c)).According to our mc-QHA calculation, the proton configuration is "more ordered" rather than "more disordered" in this pressure range corresponding to 11-18 GPa QHA .The absence of anharmonic and tunneling effects in our calculations and the problematic DFT description of H-bonds seem to constrain unrealistically the pressure stability field of HOC-11* up to 11.5 GPa EXP .Experimental observations beyond this pressure correlate at best qualitatively with our mc-QHA results at higher nominal pressures.
In the ∼11.5-18GPa EXP pressure range, d(O•••H) and d(O-H) bond lengths behave anomalously (see Fig. 6).Their pressure dependences are abruptly reversed after reaching limiting values at ∼11.5 GPa EXP .However, they resume their normal pressure-induced behavior, i.e., decrease in d(O•••H) and increase in d(O-H) with pressure soon after.OH-bending mode frequencies also have different compressive behavior below and above ∼11.5 GPa EXP [14] (see Fig. S4(a)).These changes suggest a proton configuration change at this pressure followed by normal but accelerated bond length compressive behavior, terminating in H-bond symmetrization at 18 GPa EXP .Mixing of OH-bending and stretching modes in HOC-11* (anharmonicity) starts at ∼18 GPa QHA (∼11.5 GPa EXP ) (see Sec. III B, Figs.3(b,f,j)), indicating this latter anomalous stage involves anharmonic effects and is likely accompanied by tunneling.Although tunneling is anticipated close to H-bond symmetrization, it is also expected near configuration changes.HOC-11* becomes unstable at ∼20 GPa QHA (∼14 GPa EXP ), and only HOC-22 and -11 configurations, both with fully aligned protons and with the same multiplicity (see Fig. 1) survive after that (see Fig. 4(a)).Both configurations have slightly larger volumes than HOC-11* and might explain the decrease in compressibility in this pressure range (see Figs. 4(b,c)).Suppose these configurations are present in this pressure range.In that case, mc-QHA results predict the reappearance of the highest OH-stretching frequency modes or rather an increase of the OH-stretching band frequency (compare Fig. 3(b) with Figs.1(a,d)), which has not been reported yet.Coherent proton tunneling within some length scale between these configurations is consistent with pre-H-bond symmetrization behavior.

Figure 2 .
Figure 2. (a) Volume dependence of NEB activation energies, or energy barrier heights for a single (but periodically repeated) proton jump from HOC-1 to HOC-2.Dashed and solid lines indicate activation energies for forwarding and backward jumps; colors denote the supercell size in the NEB calculation.Green lines generally overlap with yellow ones.The bottom and top axes show the unit-cell volume (2 f.u.) and corresponding experimental pressure and static PBE pressure of the base model (HOC-1) for reference.(b) Base and jumped structures used in the NEB calculations for different supercell sizes are shown in a 2 × 2 × 4 supercell.The supercells are constructed by stacking HOC-1 (base) and HOC-2 (jumped) unit cells indicated by white and orange cubes.HOC-1 and HOC-2 configurations differ by a single change in proton configuration, i.e., a single proton jump per primitive cell.

Figure 3 .
Figure 3. Pressure evolution of phonon properties.(a-d) H-projected vibrational density of states (H-VDoS) vs. unit-cell volume.The transparency of color indicates VDoS amplitude; (e-h) Γ point phonon frequencies of OH stretching (red), in-plane bending (blue), and out-of-plane bending (green) modes vs. unit cell volume.The transparency of color indicates projected strength; (i-l)angle between H-displacement vector and OH bond direction vs. volume.Symbols are from PBE results, and lines are spline interpolations to help identify mode-continuity with volume.The main bottom axis indicates the unit cell volume (2 f.u.); Extra axes at the top show the experimental and PBE static pressure at the investigated volumes.The additional axis at the bottom, d(OO) from static PBE calculations, is given to facilitate discussions.Scattered "+" and "×" symbols in panels (a-h) denote the frequency of IR and Raman peaks extracted from Refs.[27] and[13].

Figure 5 .
Figure 5. (a-d) Diffraction pattern at d hkl = 0.6-3.5 Å for 4 HOC supercell structures at 0 GPa PBE .The diffraction intensity is normalized to 1.The shaded panel is a zoom-in view at the 021 peak region (d hkl = 1.6-1.8Å, intensity amplified 10x).The arrow points to the expected d hkl position (∼1.7 Å) of the 021 peak.(e) Relative diffraction intensities of 021 peaks normalized by that of the 110 peak for HOC-11, -22 (colored curves), and "overall" (bold black solid curves) vs. volume compared.Open circles are data reported in Ref.[16].The relative 021 peak intensities for HOC-11* and -22 are less than 0.002 and, therefore, not plotted.The "overall" peak intensity is the average of peak intensities from different configurations using the mc-QHA populations at 300 K. Diffraction patterns, peak positions, and peak intensities were calculated using CrystalDiffract®[53] and the optimized crystal structures.

Figure 6 .
Figure 6.Interatomic distances d(OO) and d(OH) vs. unit-cell volume (2 f.u.) for four δ-AlOOH configurations indicated by color.d(OH)s larger and smaller than 1.2 Å correspond to H-bonds and ionic bonds, respectively.Black lines denote d(OH)'s averaged by population (nm).Scattered symbols denote d(OO)'s and d(OH)'s reported in Ref.[16].Additional axes showing 300 K experimental pressure, 300 K mc-QHA pressure, and the experimental d(OO) at corresponding volumes are also offered to facilitate comparison.Volume in the bottom axis corresponds to 2 f.u..

Table I .
Summary of 300 K phase changes in δ ] • OH-bending and stretching mode mixing (anharmonicity) likely accompanied by tunneling in HOC-11* starting at 11.5 GPa EXP (or ∼18 GPa QHA ) • HOC-11* becomes vibrationally unstable at ∼20 GPa QHA (or ∼14 GPa EXP ) • Unstable mode in HOC-11* expected to have a tunneling component near 20 GPa QHA • Predicted increase of OH-stretching bond frequencies 18 GPa EXP Symmetric H-bond • Single OH bond length • Frequency vs. volume slope changes again in OH-bending modes [27]