Probing the state of hydrogen in δ -AlOOH at mantle conditions with machine learning potential

,


I. INTRODUCTION
H-bearing mineral phases are responsible for water circulation from the Earth's surface to the interior via subduction and convection.Up to ten oceans could be stored in the Earth's interior as hydrous phases or nominally anhydrous minerals (NAMs), where H exists as substitutional or interstitial defects [1,2].Hydrogen bonds (H-bonds) transform under elevated pressures and temperatures in these solids.Like H-bonds in H 2 O-ice, they symmetrize under pressure, producing ionic bonds [3,4] or disorder before melting [5,6].Superionic diffusion of protons is also expected as a precursor to dehydration reactions responsible for melt production and volcanism.In NAMs, the H-concentration changes phase relations and plastic properties, causing seismological properties irregularities in the mantle (e.g., [7][8][9][10]).Fragile H-bonds and hydrous defects weaken the rock's rheological properties, facilitating plastic deformation and thermal convection in the mantle, a central process in Earth's evolution (e.g., [11][12][13]).Therefore, a detailed understanding of H-bond behavior and related mineral properties is essential for understanding the Earth's interior's geochemical activity, dynamic behavior, and seismological properties.
Because the pressures and temperatures (P, T s) in the Earth's interior challenge the experimental determination of materials' properties, insights from ab initio studies have been indispensable.Standard LDA [14] and PBE [15] functionals combined with phonon calculations and the quasiharmonic approximation (QHA) [16][17][18] or with ab initio molecular dynamics (MD) [19] have been used to address the physical properties of such phases.However, anharmonicity, small MD simulation size, inadequate DFT functionals, quantum nuclear effects, and other factors prevented the reliable description of H-bond behavior and property changes in these phases at relevant conditions.
Here, we focus on δ-AlOOH (δ) [20], a prototypical hydrous phase stable throughout the entire pressure range of the Earth's mantle (up to 135 GPa) [21].It has been extensively studied experimentally (e.g., [21][22][23]) and computationally (e.g., [24][25][26][27]).Therefore, experimental uncertainties and ab initio calculations' limitations to reproducing observations are well established.This phase has a simple high-pressure CaCl 2type structure [20,28] and exhibits anomalous elastic properties under pressure [26], with proton diffusion expected at high temperatures [29].δ, its isostructural siblings ϵ-FeOOH, MgSiO 4 H 2 phase H [30], the NAM CaCl 2 -type SiO 2 , along with their multiple solid solutions, are the main H-bearing phases in the deep mantle [31].Although δ-AlOOH is not the most dominant hydrous phase in the lower mantle or core, the structural simplicity, the abundance of careful measurements and calculations of δ's properties, and its prototypical nature make it a most suitable model phase for testing the performance of ab initio-based machine learning (ML) methods in addressing these systems at high P, T s.
According to measurements, H-bonds in δ symmetrize and form H-centered (HC) bonds at ∼18 GPa [22,23].Quasiharmonic approximation (QHA) based methods (e.g., [16][17][18]), while widely used to describe finite-temperature properties of solids, rely on stable phonon modes.While these methods describe effects caused by the anomalous pressure dependence of OH-stretching modes observed in the infrared (IR) spectroscopy below 12 GPa [27,32], they cannot describe δ's properties in the pressure range of H-bond symmetrization, a transition accompanied by strong anharmonicity and an arXiv:2309.06712v4[physics.comp-ph]12 Mar 2024 order-disorder precursor below ∼40 GPa in static calculations [25,27].Besides, previous ab initio QHA studies have used the LDA and PBE/GGA functionals, making comparing ab initio predictions and experimental measurements at the same pressures challenging.The relatively low experimental H-bond symmetrization pressure of ∼8-18 GPa has not been correctly reproduced in these DFT calculations because of the process's complexity and these functional's inadequate description of Hbonds, giving a transition pressure in the 30-40 GPa pressure range [23,27].
At high temperatures, superionic proton diffusion is anticipated in H-bearing systems (e.g., antigorite, diopside with H defects), usually as a precursor to dehydration.Diffusion, critical to understanding dehydration dynamics, has also been found in other H-bearing systems like antigorite [33], diopside with H defects [34], and FeOOH [35] and has been extensively investigated in ice before melting (e.g., [36][37][38][39][40][41]). A recent Born-Oppenheimer molecular dynamics (BOMD) study using the PBE functional reported proton diffusion in δ at 2,700-3,000 K [29], a temperature range higher than the 1,600-2,500 K dissociation (2 AlOOH → H 2 O + Al 2 O 3 ) temperature measured in the 20-140 GPa range [21,28,42].The PBE functional and the small MD simulation cell sizes with only a few hundred atoms are likely responsible for this discrepancy.In summary, accurately reproducing P, T stability fields of hydrous phases and their properties is fundamental for mapping out phase relations and dehydration sequences along a specific geotherm (e.g., [11]).These results provide the basis for interpreting irregular seismic properties caused by the non-homogeneous distribution of these H-bearing solids and melts in the mantle and estimating the water content in Earth's interior (e.g., [1,2,12]).
This challenge can be overcome by adopting ML-based potentials (e.g., [43][44][45][46][47][48]) in MD simulations with DFT predictive power but significantly reduced computational cost [49].This approach only invokes DFT while training the potential to reproduce interatomic forces and energies using ML-based descriptors.The ML-based MD simulations are then performed at a cost and scaling close to empirical force-field simulations.Active learning schemes (e.g., [50,51]) help create compact reference datasets by iteratively discovering and actively adding "missing" necessary atomic configurations, reducing the DFT computation cost in the potential training process.In particular, Deep Potential (DP) is a neural-network (NN) type ML potential [43] with a neural network's flexible fitting capability, allowing NN-potentials to represent chemical systems of varied nature.Recent benchmarks have shown that DP forces and energies in solids and liquids structures trained on a few thousand reference configurations can be highly accurate [49,52,53].The adoption of ML-based potentials allows the use of the strongly constrained and appropriately normed (SCAN) metageneralized gradient approximation (meta-GGA) functional's more accurate description of the H-bonded systems [54] regardless of being computationally more expensive.This method successfully calculated the water/ice phase diagram [55] and proton diffusion in liquid water [56].
Combining deep learning potentials with the SCAN functional is a promising path to simulate H-bearing systems at extreme geophysical conditions accurately.Here, we apply this hybrid approach to δ, a prototypical H-bearing system, to understand how well it overcomes the limitations presented in conventional purely DFT studies.This method enables us to perform long and large simulations that densely cover a wide (P, T ) or (V, T ) range with changing H-bond or proton diffusion behavior.We benchmark the potential against SCAN calculations and experimental measurements to assess the accuracy level achieved.We address δ's compression curve in its high-P H-bond symmetrization regime and its high-T proton diffusion regime, which were not accurately described in previous QHAor BOMD-based studies.FIG. 1. LDA, PBE, and SCAN static compression curves compared to 300 K measurements [21,22].
We employ the SCAN functional [54] to describe δ's potential energy (Born-Oppenheimer) surface and interatomic forces.Compared to LDA's underestimation and PBE's overestimation of the pressure, SCAN's static compression curve is the closest to the 300 K measured one (see Fig. 1).The slight underestimation of the SCAN volume is eliminated when thermal effects at 300 K are included, as shown in the subsequent discussion.Quantum nuclear effects disregarded in the present study might spoil this good agreement [57], and this effect must be further inspected.We are interested in the temperature range of typical subducting tectonic plates (slabs), i.e., along "slab geotherms" [58,59], where this phase is stable in the mantle, i.e., ∼ 1,000 K < T < 3,000 K [21,60].At these temperatures, classical MD should adequately describe the ionic dynamics in δ.
Our SCAN-DP interatomic potential trained on just a few thousand reference configurations can reproduce SCAN-DFT's forces and energies for configurations sampled from SCAN-BOMD simulations with 128 atoms, as illustrated in Fig. 2. The root-mean-square error (RMSE) [61] of the SCAN-DP predictions for these tests at various (V, T ) states are listed in Table SI.Configurations generated at higher temperatures produce larger RMSEs in general.At 3,000 K, the RMSEs for potential energy and force predictions are ∼2 meV/atom, and ∼0.12 eV/Å, respectively.Such accuracy is similar to previous benchmarks (e.g., [49]) in similar DP studies.SCAN-DPMD also reproduces SCAN-BOMD's structure and bonding properties at high-P, T , as illustrated by the comparison of SCAN-BOMD's and SCAN-DPMD's radial distribu-tion functions, g(r), in Fig. 3.The results are for 54 Å 3 volume, corresponding to 9.5-28.3GPa at 600-3,000 K obtained in 128-atom simulations.Plots detailing the contribution from each type of atomic pair at other P, T 's are shown in Fig. S1 [62].
The g(r) at V = 54 Å 3 shown in Fig. 3 is particularly significant because this volume approximately corresponds to the critical volume of the experimentally observed H-bond disorder transition at 300 K [22].At 600 K, the first peak at ∼1.1 Å corresponds to the ionic OH bond length.This peak is asymmetric at this volume and at all sampled temperatures due to asymmetric proton motion.At 600 K, g(r) shows a near double-peak distribution of OH ionic bonds and H-bonds; the broad shoulder centered at ∼1.4 Å corresponds to the distribution of H-bond lengths.The increased overlap of the two peaks with increasing temperature suggests the onset of a disordered state, a precursor to H-bond symmetrization [23].A similar OH bond-length distribution was observed in 300 K classical BOMD simulations with a quantum thermostat [63].Some features on g(r) disappear at elevated temperatures.For example, the 2.4 Å peak, a combination of Al-O and O-H contributions, and the 2.7 Å peak associated with O-O bonds (see Fig. S1) are split at 600 K but merge above 1,200 K.This phenomenon results from increased atomic vibration amplitudes, proton dynamic disorder, and diffusion onset at elevated temperatures.We will elaborate further on the diffusion process in Sec.II C.Under pressure, the first peak at ∼1.1 Å and the shoulder at 1.4 Å evolve into a single symmetric peak due to H-bond disorder and finally symmetrization [23,27] (see Fig. S1).
SCAN-DPMD with ab initio-level accuracy and improved efficiency will help predict the thermoelastic properties of hydrous solids, a property of first-order importance in geophysics [18,53,64].Here, we compare its predictions for δ's high-temperature compressive behavior with ab initio SCAN-BOMD predictions and measurements in the 300-3,000 K and 20-115 GPa range.Fig. 4(a) shows excellent agreement between SCAN-DPMD's high-temperature compression curves fitted to third-order Birch-Murnaghan equations of state (EoS) and the SCAN-BOMD's predicted P, V data points.Comparison with measurements [21,22] in Fig. 4(b) shows a promising 3 GPa maximum error in pressure prediction at the same V, T s.The next section will present a more detailed analysis of the 300 K compression curve at lower pressure and compare it with measurements.
These good agreements between the SCAN-DFT and SCAN-DP predictions of forces, potential energies, g(r), and high-P, T EoSs suggest that the deep-learning potentials have reached satisfactory accuracy and are predictive at least in the 20-120 GPa pressure range and up to 3,000 K.

B. The 300 K compression curve and dynamic stabilization of the HC phase
The 300 K compression curve is one of δ's best-determined properties.Experimentally, δ exhibits an anomalous change in compressive behavior at ∼10 GPa, which is usually attributed to a change in H-bond structure, most often with H-bond symmetrization [24,[67][68][69], but also with H-bond disorder [22,23,70].However, it has been impossible to reproduce this compressive behavior using QHA-based DFT calculations.This is because (a) static ab initio calculations without multiconfiguration analysis do not account for H-bond disorder; (b) the HC phase atomic configuration is dynamically stable above ∼35 GPa only in static GGA/PBE calculations [25].Strong anharmonicity produces unstable phonon modes at lower pressures, hindering the application of the QHA for high-temperature calculations [25,27].Using SCAN-DPMD, we optimize and equilibrate the cell shape at several pressures to obtain the static and the 300 K compression curves shown in Figs.5(a,b).With these, we can analyze the effect of atomic motion on the H-bond state, its effect on the volume, and compare the latter with measurements [21,22,65].As shown in Fig. 5(b), extrapolations of a Birch-Murnaghan EoS fit to results below and above ∼10 GPa produce divergent curves (dotted blue and red dotted lines), indicating a change in the OH-bond state at ∼10 GPa.
The 300 K isothermal bulk modulus, K T = V (∂P/∂V) T , shown in Fig. 5(c), amplifies the subtle change in compressibility across this critical point.Across this state change in the 0-20 GPa range, the OH bond-length distribution, g OH (r), evolves from a double peak to a single modal distribution, though not symmetric (Fig. 5(d)).In this pressure range, the H-bond compresses quickly while the ionic OH bond stretches, a well-known behavior of these bonds (e.g., [57,71]).δ's compressive behavior obtained from SCAN-DP calculation confirms previous PBE/GGA results [26] that this change in compressive behavior occurs in the 30-40 GPa range in static calculations (dashed line in Fig. 5(c)).
SCAN-DPMD results at 300 K account for the thermal expansion and reproduce well the 300 K experimental compression curve [21,22,65]: in the 0-15 GPa range, the difference between results and measurements is less than the differences between measurements; at higher pressures, the predicted SCAN-DPMD volume is slightly overestimated, with an error smaller than 0.25 Å 3 /f.u.Although predicted [24], it is striking to see that the inclusion of classical ionic motion lowers the transition pressure by ∼27 GPa, giving a transition pressure that agrees quite well with the experimentally measured one at ∼9 GPa (e.g., [22,69]).This 300 K result is surprising and puzzling.At this low temperature, quantum proton motion, particularly tunneling [72], should impact the H-bond behavior (e.g., [73]).
Above 10 GPa, V and K T display a smooth monotonic dependence on pressure.Table I summarizes the Birch-Murnaghan EoS parameters resulting from the fitting of measured and calculated 300 K P, V data above 10 GPa.Fitting the SCAN-DPMD results within the 10-64 GPa pressure range with K ′ 0 = 4 gives V 0 = 55.5 Å 3 , K 0 = 225 ± 2 GPa (red dotted line in Fig. 5(b)).The predicted V 0 is in excellent agreement with the measured V 0 [22] fit to the same Birch-Murnaghan EoS, 55.47 Å 3 .The predicted bulk modulus differs by ∼3% from the experimental one, 219 GPa [22], an impressive agreement for a result extrapolated to 0 GPa.
Our SCAN-DPMD results are significantly more accurate than previous PBE-BOMD results, which reported a V 0 = 58.5 Å 3 and K 0 = 183.4GPa using the Vinet EoS fit to results above 10 GPa [63].This improvement can be attributed to (a) SCAN's more accurate description of the P-V relation, (b) the denser sampling of (P, V) states over a broader pressure range, and (c) larger, longer, and better-converged simulation runs enabled by DPMD's performance leap compared to BOMD.Nevertheless, this level of agreement between classical MD results and measurements at 300 K is unexpected.
This good agreement between SCAN-DPMD results and experimental data above 10 GPa suggests that our simulations describe well the H-bond state above this pressure.As indicated in Fig. 5(d), not even at 20 GPa, the first two peaks in the g OH (r) pair distribution function have merged, indicating that H-bond symmetrization has not yet been achieved in the simulations.Therefore, the change in compressive behavior at ∼10 GPa cannot be attributed to H-bond symmetrization in this classical simulation.Fig. 6 shows in detail the evolution of g OH (r) and g HH (r) pair distribution functions from 0 to 10 GPa.Both show two clear peaks at 0 GPa and tend to merge with increasing pressure.However, only the g HH (r) peaks are fully merged at 10 GPa.g OH (r) still displays two superposing broad peaks.This indicates that H-bonds are not symmetric in these classical simulations at 10 GPa (or 20 GPa).The change in the H-bonds' state at 10 GPa seems related to a change in H-ordering, likely into a disordered state.Our previous multi-configuration PBE-QHA study [27] suggested several short-to medium-range ordered H-bonded configurations remain dynamically stable beyond the pressure where the compressive behavior changes (∼10 GPa in experiments, ∼12 GPa in PBE-QHA calculations).This result supported the notion of a more disordered H-bond state.The current results showing an asymmetric first peak in g OH (r) and the fully merged broad peak in g HH (r), and the previously identified multi-configuration equilibrium state [27] suggest that the change in compressive behavior at ∼10 GPa is associated with the emergence of a disordered state in the simulations.Static or dynamic, disorder involves proton hopping and is a precursor to proton diffusion at higher temperatures.This 300 K behavior might change significantly if the proton dynamics is addressed quantum mechanically using path integrals [73].

C. Proton diffusion at high temperatures
The DP potential allows us to systematically perform largescale DPMD simulations to investigate proton diffusion in δ over a wide P, T range.We perform 96 NVT DPMD simulations covering the pressure range of 35-140 GPa and temperature range of 1,500-3,000 K (see Fig. S2).They contain 8,192 atoms and run for 2 ns.We do not observe diffusion, structural change, or melting of the Al and O sub-lattices in this P, T range.
Proton diffusion is quantitatively characterized by protons' mean-square displacement (MSD).Protons' MSD over a simulation run time t, L 2 H (t), is defined by [74,75] where N H is the total number of protons, R H,i (τ) denotes the position of i-th proton at moment τ, and "⟨ • ⟩" denotes the ensemble average over the start time τ.A linear (non-linear) dependence of the MSD on the simulation run time, t, reflects a continuous (irregular) diffusive behavior.The self-diffusion coefficient, D H , or the diffusivity, is related to the slope of the Figs. 7(a,b) show the protons' MSD over a time span of 1 ns for V = 47 Å 3 and T = 1,500 K (60.9 GPa) and 2,400 K (67.5 GPa), respectively.MSD at several other V, T conditions are available in Fig. S3.Increasing the temperature significantly increases protons' diffusivity (see Fig. S3).Diffusion is not observed at 1,500 K.It starts between 1,800 K and 2,100 K and becomes steady at 2,400 K. Figs.7(c,d) show the corresponding proton trajectories.Fig. 7(b), shows that the protons' MSD along the c direction is approximately twice as large as those along the a and b axes, indicating they move more freely through the interstitial channels along the c direction (see Figs. 3(a,b)).A similar phenomenon has also been reported in the NAM phase hydrous Al-bearing stishovite, which has a similar CaCl 2 -type Si-O framework [76].
At 1,800-3,000 K, the relationship between log D H and 1/T is linear at all volumes and can be fitted using the Arrhenius equation (see Fig. S4), where k B is the Boltzmann constant and E a denotes the activation energy.This behavior indicates that proton diffusion in δ is a common thermally activated process.The volume dependence of the activation energy, E a , is nearly linear (Fig. S5).
Using these relationships, we interpolate D H vs. V and T using the finite temperature EoS displayed in Fig. 4. The magnitude of protons' diffusivity is shown as the background color in Fig. 8.For diffusivity corresponding to MSD > 1 Å 2 /ns, (e.g., Fig. 7), diffusion is steady, stable, and characterized by a Ref. [29] and dissociation boundaries from Refs.[21,42] are also shown.Adiabatic mantle geotherm [77] and slab geotherm [58].
linear dependence of the MSD vs. time (reddish background area); for diffusivity of < 0.1 Å 2 /ns (bluish background area), diffusion can be intermittent with non-linear dependence of the MSD vs. time and the system is considered to be solid [78].Depending on the pressure, for 1,800 K < T < 2,100 K, D H approaches 1 Å 2 /ns and starts deviating from the high-T Arrhenius behavior (see Fig. S4).This region of the diffusivity diagram is transiting from a blueish to a whitish background.At these P, T conditions, one can recognize in Fig. S3 a change in the MSD's behavior from non-linear to linear.We designate this regime as the onset of the fully superionic behavior in the P, T phase space and represent this D H = 1 Å 2 /ns diffusion boundary in Fig. 8 with a solid black curve.
The experimental dehydration boundary reported in Ref. [21] resembles our diffusion boundary.Similar measurements in Ref. [42] pinpointed two δ dehydration temperatures in two distinct pressure ranges: 1,850-1,900 K at 28-68 GPa and ∼1,959 K at 100-110 GPa.Two divergent boundaries were reported in the mid-lower mantle (68-100 GPa) due to uncertainties related to unidentified X-ray diffraction peaks and uncertain experimental conditions (Fig. 8).Their lower boundary directly connects the low-pressure and high-pressure boundaries and runs parallel to our boundary with a 200 K downward shift.It has been argued [42] that the disagreement between these studies could be due to uncertainties in granularity or the preferred orientation of samples in the high-pressure cell, which could affect the dehydration product detection sensitivity.The proximity of our diffusion boundary and the experimental dehydration boundary, particularly their similar pressure gradients, suggests that continuous proton diffusion controls the dehydration process as in antigorite [33] or diopside with H defects [34].
The diffusion boundary previously obtained with PBE-BOMD [29] largely overestimates the superionic transition temperature.It is known that supercell size could affect the superionic diffusivity calculations, often leading to overestimation of the transition temperature [79,80].The short simulation time and the PBE's inaccurate description of the H-bond could also contribute to the overestimation of the diffusion boundary.
It is also possible to investigate the change in protons' diffusivity by inspecting the constant-volume heat capacity, C V , through the P, T range shown in Fig. S6.The constant-volume heat capacity, C V , can be calculated from MD ensemble averages [81] as where k B denotes the Boltzmann's constant, N denotes the number of atoms, ⟨E 2 ⟩ and ⟨E⟩ 2 denotes the run average of energy E and its square E 2 .C V computed using this method does not exhibit a size effect for simulation cell sizes varying from 128-27,648 atoms in both the solid and diffusive states (Fig. S6).C V of an approximately harmonic solid is expected to plateau at the Dulong-Petit limit, 3NK B .Deviation from this value suggests strongly anharmonic behavior, in the present case, entering the diffusive regime.Fig. 9 shows δ's C V vs. P, T as background color.For T < 1,500 K (blueish background), C V generally follows the Dulong-Petit limit with at most 5% excess.For T > 1,800 K, C V increases quickly.Depending on the pressure, the background color turns from blue to white then red between 1,800 K and 2,100 K, indicating > 10% deviation from 3Nk B within this color change range.The D H = 1 Å 2 /ns diffusion boundary corresponds to C V ∼6% above the Dulong-Petit limit.Mea-surements of C V could be useful to determine the onset of δ's superionic state.

III. DISCUSSION
This study demonstrates the necessity of adopting a hybrid ab initio description aided by deep-learning potential at extreme P, T mantle conditions to address the properties of a hydrous system.The DP-GEN active learning scheme used to develop the potential is highly efficient, requiring only a few thousand massively parallelized DFT calculations on reference configurations, taking only a few days.The number of ab initio calculations needed to prepare the interatomic potential for a wide range of P, T conditions, ∼3,000 in this case, is far fewer than that of a typical BOMD run necessary for a single P, T sampling, ∼10 4 -10 5 .This efficiency allows us to use more accurate functionals and perform more accurate simulations using larger simulation cells, longer run times, and denser sampling of points in P, T phase space.Yet, it is still desirable to develop these potentials further to perform under reactive conditions and detect the dehydration process in a simulation.
The current DP-SCAN approach combined with classical MD reproduces measurements of δ's room temperature compression curve surprisingly well.It paves the way for adopting the path-integral approach to quantum ionic dynamics [82,83] in low-temperature simulations of δ.Zero-point motion effects on the EoS of solids at low temperatures, particularly in H 2 O-ice (e.g., [57]), are well-known.On δ, the presence of the quantum nuclear effects, e.g., tunneling, is an ongoing debate; evidence supporting its presence [3,27] and absence [84] has both been published.Such effects could affect H-bond disordering, symmetrization, and details of P, T phase diagrams of hydrous phases (e.g., [85]) but have yet to be fully explored in hydrous minerals.At typical mantle temperatures of thousands of Kelvin, classical MD combined with the SCAN functional seems to reproduce closely the thermal EoS of a prototypical lower mantle hydrous phase changing the H-bonds' state.
At mantle and subducting slab conditions (Fig. 8), it remains challenging to predict the dehydration process, particularly the transition from the superionic state to the dehydration point.The lack of detailed experimental data in this critical range of conditions aggravates understanding of this process.Electrical conductivity measurements commonly identify diffusion.Electrical conductivity comprises ionic and electronic components: ionic conductivity, according to the Nernst-Einstein relation, will be proportional to free proton concentration, which varies exponentially with temperature; electronic conductivity would require more complex electronic structure calculation, which extends beyond the scope of this study.δ's electrical conductivity has only been measured up to 1,200 K and up to 20 GPa, and has shown exponential dependence vs. reciprocal temperature [86] similar to ours.Extending such measurements to more extreme conditions is desirable to understand better the relation between the superionic state and the dehydration process.
Simulations of these processes will be fundamental to understanding water circulation in Earth's interior.Subducting slabs carry hydrous phases into the mantle.Describing and predict-ing the behavior of such phases at extreme conditions is key to understanding mantle dynamics, e.g., volcanism, melt generation, etc.While here we focus on the superionic state, we see a relationship between the onset of the superionic state and the dehydration boundary.Our predicted "diffusion boundary" (see Fig. 8), the dehydration boundary reported by Ref. [21], and Ref. [42]'s lower boundary, all share a similar slope that is less steep than the temperature profiles along the subducting slab.These diffusion and dissociation boundaries intercept the subducting slab geotherm [58,59] at a depth of ∼2,400-2,700 km near the bottom of the lower mantle, and the normal mantle geotherm [77] at ∼1,200-1,500 km, i.e., approximately mid-mantle.These results confirm previous suggestions that δ could remain stable up to near the bottom of the mantle and only then release water.H 2 O released from δ should react with its environment to form other H-bearing phases or melts, e.g., ϵ-FeOOH and FeH x [42,87] in the deep mantle.Understanding the entire diffusion to dehydration process and the P, T conditions for these transitions will help clarify the deep water cycle, the potential signature of water presence in seismic tomography, and the effects of water on mantle properties.

METHOD Machine learning potentials
The NN potential for δ was developed based on the Deep Potential Smooth Edition (DeepPot-SE) model [48] implemented in DeePMD-kit v2.1 [88,89].Two-body embedding with coordinates of the neighboring atoms (se e2 a) was used for the descriptor.The embedding network shape is (25,50,100).The fitting network shape is (128,128,128).The cut-off radius is 6 Å, and the smoothing parameter is 0.5 Å.
The model was trained using the Adam optimizer [90] for 1× 10 6 training steps, with the learning rate exponential decaying from 1 × 10 −3 to 3.51 × 10 −8 throughout the training process.The loss function L(p e , p f ) is [88] where p e decays linearly from 1.00 to 0.02, and p f increases linearly from 1 × 10 0 to 1 × 10 3 throughout the training process.

Active learning scheme
The DP-GEN concurrent learning scheme [91] was employed to create the training data set and to generate the potential.We randomly extracted 118 labeled configurations from 59 BOMD runs at various V, T s to generate the initial potentials and to kickstart the DP-GEN training process.6 DP-GEN iterations were performed to explore the configuration space and to eventually generate a potential reaching satisfactory accuracy requirement for DPMD for a temperature range of 300 < T < 3, 000 K. for a few thousand timesteps.After the simulations, the error estimator (model deviation), ϵ t , is calculated every 50 MD steps based on the force disagreement between the candidate DPs [51,91]: where F w,i (R t ) denotes the force on the i-th atom predicted by the w-th potential for the R t configuration.For a particular configuration, if ϵ t satisfies ϵ min ≤ ϵ t ≤ ϵ max , the corresponding configuration is collected, then labeled with DFT forces and total energy then added into the training dataset; if ϵ t < ϵ min , these configurations are considered "covered" by the current training dataset; if ϵ t > ϵ max are considered failed and were discarded.After a few iterations, almost no new configurations are collected according to this standard (> 99% are "accurate" for a few iterations), the DP-GEN process is then complete.The parameters for these DP-GEN iterations are listed in Table SII in detail.After these DP-GEN iterations, our training dataset consists of 3,487 configurations labeled with ab initio force and energy.DPMD simulations were performed using the LAMMPS code [92] with 0.5 fs timesteps.

DFT calculations
Training and testing data sets were based on ab initio calculations performed with the Vienna Ab initio Simulation Package (VASP) v6.3 [93].The strongly constrained and appropriately normed (SCAN) [94] meta-GGA functional with PAW basis sets were adopted.The cutoff energy for the plane-wave-basis set was set to 520 eV.The Brillouin zone sampling for the 2 × 2 × 4 supercells (with 16 f.u. of δ, or 128 atoms) used a shifted 2 × 2 × 2 Monkhorst-Pack k-point mesh.BOMD simulations were performed with a timestep of 0.5 fs.
The investigation of the 300 K compression curve involves DPMD simulations with 1,024-atom (or 4 × 4 × 8) supercells.We perform NPT simulations for 50 ps to equilibrate the structure.Then, perform NVT simulations for another 50 ps to obtain the pressure at the given volume.The investigation of high P, T s, involved systematic DPMD simulations with 8,192atom (or 8 × 8 × 16) supercells.We perform NVT simulations.96 DPMD simulations at different (P, T )'s cover the P, T range from 1,500 K to 3,000 K and 10 GPa to 180 GPa (see Fig. S2).
Each MD simulation at equilibrated V, T conditions runs for 2 ns, or 4 × 10 6 timesteps.These simulations were performed concurrently on 96 GPUs.

Workflow management
Our workflows, including the DP-GEN active learning iterations and subsequent MD simulations, were implemented and managed using the Snakemake [97,98]

FIG. 7 .
FIG. 7. Proton diffusion at high temperatures for V = 47 Å 3 conventional cell which corresponds to 60.9 and 67.5 GPa at 1,500 K and 2,400 K. (a) Proton mean-square displacement (MSD) at 1,500 K and (b) 2,400 K. (c) Proton trajectories at 1,500 K and (d) 2,400 K; blue, red, and gray dots denote Al, O, and H ions.
FIG. 8.The diffusion phase diagram of δ.The background color represents the diffusivity D H .The solid black curve indicates the D H = 1 Å 2 /ns boundary for steady proton diffusion characterized by a linear dependence of MSD vs. time.The diffusion boundary fromRef.[29]and dissociation boundaries from Refs.[21,42] are also shown.Adiabatic mantle geotherm[77] and slab geotherm[58].

B
FIG. 9.The ratio between C V and the Dulong-Petit limit, 3Nk B , at different pressures and temperatures.Black symbols indicate a 6% excess above the limit.The error bar is determined based on the P, T grid size.Results are compared to the D H = 1 boundary and measured dehydration phase boundaries[21,42].

FIGURE S6 -
FIGURE S6-Size effect on C V at ∼110 GPa.

TABLE I .
Comparision of 300 K EoS parameters for HC-δ.
workflow management system.

TABLE SI -
Details for RMSE of force and potential energies for validation.

TABLE SII -
Details of the DP-GEN iterations for generating the training dataset.