Edinburgh Research Explorer Plastic and Superionic Helium Ammonia Compounds under High Pressure and High Temperature

Both helium and ammonia are main components of icy giant planets. While ammonia is very reactive, helium is the most inert element in the universe. It is of great interest whether ammonia and helium can react with each other under planetary conditions, and if so, what kinds of structures and states of matter can form. Here, using crystal structure prediction methods and first-principles calculations, we report three new stable stoichiometries and eight new stable phases of He-NH 3 compounds under pressures up to 500 GPa. These structures may exhibit perovskitelike structures for HeNH 3 and He 2 NH 3 , and a host-guest crystal structure for He ð NH 3 Þ 2 . Superionic states are found in all these He-NH 3 compounds under high pressures and temperatures in which the hydrogen atoms are diffusive while the nitrogen and helium atoms remain fixed. Such dynamical behavior in helium ammonia compounds is quite different from that in helium water compounds, where weakly interacting helium is more diffusive than stronger bound hydrogen. The low-density host-guest phase of space group I 4 cm is found to be stable at very low pressures (about 3 GPa) and it enters into a plastic state, characterized by freely rotating ammonia molecules. The present results suggest that plastic or superionic helium ammonia compounds may exist under planetary conditions, and helium contributes crucially to the exotic physics and chemistry observed under extreme conditions.


I. INTRODUCTION
Helium is the second most abundant element in the universe, and it occupies a large fraction of the atmospheres of Uranus and Neptune. Although helium is generally considered to be unreactive, several groups have reported stable helium compounds under high pressure. For instance, it can react with metals [1,2], metal oxides, sulfides, or fluorides [3][4][5], in which the insertion of helium atoms significantly modifies the long-range Coulomb interactions [5]. In addition, helium can react with inert gases [6][7][8], nitrogen [9,10], water [11][12][13], and ammonia [14], where it exhibits negligible charge transfer and plays a role in increasing the internal pressure.
Since the pioneering work of predecessors [15][16][17], considerable effort has been devoted to studying superionic states in hot dense molecular compounds, which are characterized by almost free diffusion of protons through fixed oxygen or nitrogen sublattices. Several theoretical [18][19][20][21][22][23][24][25] and experimental studies [26][27][28][29][30] have investigated the constituents of celestial bodies, which has led to theoretical models for planetary interiors that may exist within the "hot ice" layers of Uranus and Neptune. The hot ice is composed of water, ammonia, and methane under extreme pressures up to several megabars and temperatures up to several thousand kelvin, and is below the atmospheres of hydrogen and helium. Systematic computational investigations of these mixtures have achieved significant advances [31][32][33][34][35][36]. For instance, ammonia was predicted to react with water and form a series of compounds (with compositions ranging from 4∶1 to 1∶2) along the isentropes of Uranus and Neptune [33]. Bethkenhagen et al. [37] studied the 1∶4∶2 ternary mixture of ammonia, water, and methane, which is close to their solar abundance ratios, and constructed a new adiabatic model of Uranus. These studies of the mixtures in the hot ice layer have shed light on the understanding of these icy planets. This layer is likely to interact with the constituents of the upper planetary atmospheres; therefore, it is essential to understand how the molecular ices react with hydrogen or helium at high-pressure and high-temperature conditions. For instance, a recent study reported unusual chemical reactions between ammonia and hydrogen under those conditions [38]. The exotic ionic mobility behavior seen in superionic states is mirrored in applications in material science, where it greatly enhances the performance of lithium-based battery materials [39][40][41] and copper-based thermoelectric materials [42,43].
Apart from superionic states in hot dense molecular ammonia, the plastic phase in ammonia, in which protons rotate around the fixed nitrogen atoms, has been reported to emerge below the superionic pressure range [27]. The plastic phase and rotational motion in water have been studied under high pressure [44][45][46]. Very recently, Li et al. [47] proposed that the plastic phases of materials with rotational disorder may be used as a new refrigeration agent, which absorbs or releases a large amount of energy when the systems enter or leave the plastic region through tuning of pressure.
Both helium and ammonia are major components of icy giant planets, but whether they can form new stable compounds under pressure and whether plastic and superionic states can exist in these compounds are still open questions. Also motivated by our previous study of the multiple superionic states in helium water compounds [13], we aim to explore the formation, structural evolution, and dynamic properties of helium ammonia compounds under high-pressure and high-temperature conditions. We have studied the helium ammonia system using crystal structure searching methods and have found eight new stable phases, some of which exhibit interesting host-guest and perovskitelike structures. We also found interesting plastic and superionic states in these compounds, which behave very differently from the helium water system.

II. RESULTS
We have explored the crystal structures in helium ammonia compounds below 500 GPa using variablecomposition structure prediction methods in conjunction with density functional theory (DFT) total energy calculations. A structure is considered to be thermodynamically stable if its enthalpy of formation is negative relative to either elements or other potential compounds. We used the most stable solid phases of the related elements or compounds at the pressure range we studied to calculate the formation enthalpy. As the convex hulls shown in Fig. 1(a) and the phase diagram shown in Fig. 1(b), we found three stable compositions and eight new stable structures in the pressure range studied. These structures belong to three different stoichiometries: HeðNH 3 Þ 2 , HeNH 3 , and He 2 NH 3 . We have checked the stability of the newly predicted phases against different exchangecorrelation functionals including different van der Waals (vdW) corrections, hybrid functionals, and the random phase approximation (RPA), as shown in Figs. S1-S3 in Supplemental Material [48], which all demonstrate that the results are qualitatively independent of functional choice. Phonon calculations indicate that the above structures are all dynamically stable, as shown in Fig. S4 of Supplemental Material [48]. The stable pressure range of the phases of each composition are shown in Fig. 1(b), and the pressure-induced phase transition sequences are as follows: Pnma ⟶ 165 GPa P2 1 2 1 2 1 ; It should be noted that there are several energetically competitive structures for these type of molecular crystals. For instance, the enthalpies of Pbcn and Cmc2 1 HeðNH 3 Þ 2 structures are very similar to those of the I4cm phase at 10 GPa. A Pca2 1 HeðNH 3 Þ 2 structure has an enthalpy close to that of the Fmm2 HeðNH 3 Þ 2 phase at 100 GPa. A Pna2 1 HeðNH 3 Þ 2 structure is only about 1 meV=atom higher than the Pbcm HeðNH 3 Þ 2 phase at 300 GPa, while a Ima2 HeðNH 3 Þ 2 structure is only about 2 meV=atom higher in enthalpy than the Cm HeðNH 3 Þ 2 phase at 400 GPa. The lattice parameters for all the competitive structures are listed in Table I in the Supplemental Material [48].
Interestingly, we found that helium ammonia compounds with helium/ammonia ratios of 1∶1 and 2∶1 are inclined to form perovskitelike structures in which each nitrogen atom resides within an octahedron formed by six hydrogen atoms, which are connected to the nitrogen atoms via three covalent bonds and three hydrogen bonds, as shown in Fig. 1(c). Because of the different bond length between the covalent and hydrogen bonds, the ammonia octahedra are highly distorted. In Pnma HeNH 3 , half of the voids are empty and the other half are doubly occupied, while the voids in P2 1 =c He 2 NH 3 are all doubly occupied due to the increased helium atom density, as can be seen in Fig. S5 of the Supplemental Material [48]. According to their enthalpies, our predicted Pnma phase HeNH 3 is more stable than the previously reported P2 1 2 1 2 1 phase [14] at pressures below 165 GPa, and it extends the stability range of HeNH 3 down to 28 GPa.
Different from perovskitelike structures, all compounds in the HeðNH 3 Þ 2 stoichiometry belong to a host-guest prototype. As shown in Figs. 1(d)-1(f), looking along one axis, all of the nitrogen atoms are arranged in four-and eight-membered rings, and the helium atoms are distributed at the center of the eight-membered rings. In particular, the I4cm phase becomes stable at quite low pressures (about 3 GPa). At higher pressures, the ammonia molecules in HeðNH 3 Þ 2 compounds are ionized by compression and form various chemical units such as NH, NH 2 , NH 3 , NH 4 , and even long N-H chains. More crystal structures of HeðNH 3 Þ 2 of this type can be found in Fig. S6 of Supplemental Material [48]. From our band structure calculations, both at PBE-D3 [49] and hybrid PBE0 level of theory [50], all of the helium ammonia compounds are insulating with wide band gaps, as shown in Figs. S7 and S8 [48]. The I4cm HeðNH 3 Þ 2 structure has a tunable band gap and an abnormal pressure dependence with a sharp increase in the band gap of around 0.064 eV=GPa over the pressure range 3-34 GPa. Recently, Liu et al. [34] carried out experimental and theoretical studies on water ammonia mixtures at extreme conditions and found that the intermediate products can be geometrically protected and thus survive for long times, which may also be a challenge for the future experimental verifications of our prediction on pure helium ammonia compounds.
To study the dynamical properties of the predicted helium ammonia compounds, we start from the HeNH 3 composition, which is stable in a relatively large and continuous pressure range. We performed extensive ab initio molecular dynamics (AIMD) simulations based on the Born-Oppenheimer approximation within the pressure range of 100-500 GPa up to the melting temperature. Diffusion coefficients (D) were calculated for the nitrogen, hydrogen, and helium atoms as calculated from the averaged mean-square displacements (MSD) to monitor different types of atomic motions. More details about the methods can be found in Sec. V. This allowed us to distinguish the solid, plastic, superionic, and fluid phases.
First, we analyzed the atomic trajectories from simulations at the initial pressure of 100 GPa for the Pnma HeNH 3 phase. We calculated the averaged MSD [Figs. 2(a) and 2(b)] and compared the atomic motions in trajectories [Figs. 2(e) and 2(f)] along one axis of the simulation box. We found that the atoms oscillate around the equilibrium positions at 500 K (D ¼ 0), while hydrogen becomes diffusive within the static nitrogen and helium frameworks . This clearly differentiates two different states: the solid phase (500 K) and the superionic phase (2000 K). Furthermore, we analyzed the atomic motions in such phases in terms of the radial distribution function (RDF), as shown in Fig. S9(a) of the Supplemental Material [48]. By heating to 2000 K, each of the RDFs related to hydrogen became more liquidlike, while the RDFs that are not related to hydrogen remained solidlike.
Inspired by the temperature-induced features in helium ammonia compounds discussed above, we expanded the pressure range to explore the superionic region up to about 500 GPa with a pressure increment of 100 GPa, as shown in Fig. 3(a). Every colored symbol corresponds to an independent simulation to avoid memory effects, and the pressures and temperatures were obtained as statistical averages. According to the averaged MSD, RDF, and the averaged structure, the final states of these simulations were classified into four types: solid (blue squares), plastic (dark green diamonds), superionic (cyan triangles), and fluid (orange circles). Dashed lines were fitted from the boundary points of each type of atomic motion which divide the phase diagram into four regions: the solid, plastic, superionic, and fluid phases. On the boundary between the superionic and fluid phases we found that helium atoms begin to become diffusive together with superionic hydrogen within a fixed nitrogen sublattice, as shown in Fig. S11 [48]. However, the diffusion coefficient of helium is almost 2 orders of magnitude lower than that of hydrogen, which is very different from superionic helium water compounds in which the helium atoms diffuse faster than the protons [13]. Since NVT simulations tend to overestimate melting temperature, such a helium diffusion regime in helium ammonia compounds may actually indicate an approach of the melting temperature; however, the possible superionic helium region should be very small if there is any. Compared with pure ammonia, the insertion of helium greatly extends the superionic region to much higher temperatures above 200 GPa; see temperature than the stability region of superionic helium ammonia compounds, it is possible that superionic helium ammonia might still form in interiors of other giant icy planets with different isentropes due to differences in mass, composition, or the distances to their main stars.
To explore the possible plastic state of the helium ammonia compounds at low pressures, we performed AIMD simulations for the host-guest I4cm HeðNH 3 Þ 2 phase. As expected, we found that the I4cm structure remains solid at low temperature, as shown in the MSD and atomic trajectory in Figs. 2(c) and 2(g). However, it enters a plastic state at around 10-40 GPa and 500-1200 K, in which hydrogen atoms rotate around the fixed nitrogen atoms, as shown in Fig. 3(b). One can find an atomic trajectory for such a plastic phase in Fig. 2(h), in which the three main distribution areas of hydrogen atoms surrounding the nitrogen atoms overlap due to the rotation of the ammonia molecules. The crucial signature of this type of plastic phase is a finite MSD with a plateau at a quite low value, as shown in Fig. 2(d). In the RDFs, the signature of hydrogen atoms in the plastic phase looks liquidlike, which is not very different from the superionic phase, as shown in Fig. S9(b) of the Supplemental Material [48].
On the boundary between the plastic and fluid phases we found in some trajectories that helium atoms become diffusive around the freely rotating ammonia molecules, as shown in Fig. 3(b) and Fig. S12 [48]. However, the pressures and temperatures of the plastic phases were lower than those of the superionic phase. Note that if we model the Pnma HeNH 3 phase at lower pressures and temperatures, it also forms a plastic state, which has been outlined in Fig. 3(a). We find that both compounds exhibit a negative slope for the phase boundary between the plastic and superionic phases in the pressure-temperature phase diagram, which also occurs in pure water ice [25]. More details can be found in Figs. S13 and S14 in the Supplemental Material [48].

III. DISCUSSION
Hermann et al. [51] and French et al. [23] found that nuclear quantum effects significantly affect the water ice phase diagram at megabar pressures. The isotope effect may also affect the dynamical properties including the diffusion coefficient and vibrational spectrum. We therefore studied this in our helium ammonia system and carried out two independent AIMD simulations for the superionic states of HeðNH 3 Þ 2 and HeðND 3 Þ 2 at around 2000 K and 100 GPa. We have estimated the nuclear quantum effects by comparing the vibrational density of states, as shown in Fig. 4. In the harmonic zero temperature phonon DOS, the high-frequency phonon peak decreases from about 3600 cm −1 for hydrogen to about 2600 cm −1 for deuterium, which is close to the ratio of 1= ffiffi ffi 2 p . As it was found that there is almost no charge transfer between helium and ammonia, such changes in the phonon frequencies are absent in the helium phonon DOS. When increasing the temperature to 2000 K, the vibrational frequencies of both nitrogen and helium are almost identical in HeðNH 3 Þ 2 and HeðND 3 Þ 2 , and they also agree with their phonon DOS at zero temperature. However, both H and D have a nonzero DOS at zero frequency, which represents a finite diffusion coefficient. Their highest frequency modes are related to a redistribution of the N─HðN─DÞ bonds length, which represents the emergence of the superionic hydrogen/deuterium state. We calculated the diffusion coefficient from the averaged MSD, which decreased with increasing mass from 1.4 × 10 −8 m 2 =s for H to 0.88 × 10 −8 m 2 =s for D. The appearance of superionic phases in both HeðNH 3 Þ 2 and HeðND 3 Þ 2 indicates that nuclear FIG. 3. Proposed phase diagram of HeNH 3 (a) and HeðNH 3 Þ 2 (b) at high pressures obtained from structure searches and AIMD simulations. The simulations are marked with four different symbols: blue square, dark green diamond, cyan triangle, and orange circle represent the solid, plastic, superionic, and fluid states, respectively. Black dashed lines are fitted to the phase transition boundaries. The calculated isentrope for Uranus (Neptune) is shown as a thick dark green (dark blue) line. The white dash-dotted line covers the superionic phase in pure ammonia, taken from Redmer et al. [20]. Red solid circles mark some trajectories in which the superionic hydrogen or plastic phases coexist with diffusive helium, which has much smaller diffusion coefficients than superionic hydrogen near the phase boundaries.
quantum effects do not affect the existence of the superionic state in the helium ammonia compounds, although we do not know whether the nuclear quantum effects affect the onset of the superionic state. We also compared the internal energy with the energy contribution from the nuclear quantum corrections, and found the influence to indeed be small, as listed in Table II in the Supplemental Material [48].
To cross-check our results of nuclear quantum effects, we have carried out path integral molecular dynamics (PIMD) with a colored noise approach, which has been successfully used in recent works by Ceriotti et al. [52,53] and Bronstein et al. [54]. We calculated the averaged meansquare displacement of the trajectory from such advanced PIMD simulation with eight beads, shown as Fig. S15 in Ref. [48]. Although the diffusion efficient is slightly reduced from 1.417 × 10 −8 m 2 =s (original AIMD) to 1.05 × 10 −8 m 2 =s (PIMD with colored noise approach), the conclusion remains the same that the superionic state of hydrogen atoms still exists in the Pnma HeNH 3 under 120 GPa and 2000 K.
Our previous work showed that helium can become freely diffusive in a fixed water ice sublattice in helium water compounds and the diffusion coefficient of helium can even be higher than that of superionic hydrogen [13]. This can be explained because the helium atoms are weakly interacting with the ice network, and hence can move freely; on the other hand, protons are bound in oxygen with covalent bonds and require more energy to become diffusive. Strangely, such a superionic helium phase is absent in the helium ammonia compounds studied here. Although on the boundary of the superionic phase and fluid phase in helium ammonia compounds, some trajectories show that the helium atoms tend to be diffusive, the diffusion coefficient of helium in those cases is much lower than that of superionic hydrogen. This is very different from the case of the helium water system. To understand the reason for this difference, we carried out reduced density gradient (RDG) analysis [55] to investigate the interactions in the Pnma phase HeNH 3 at 100 GPa, and compared the results with the I4 1 md He 2 H 2 O phase at 10 GPa and the Fd3m He 2 H 2 O phase at 100 GPa, as shown in Figs. 5(a)-5(c). The areas surrounded by solid magenta and dashed green lines represent the hydrogen bonding and vdW interactions in the different structures, which clearly show that the hydrogen bonding and vdW interactions in the Pnma phase of HeNH 3 are more significant and discrete than in the I4 1 md He 2 H 2 O phase at 10 GPa and the Fd3m He 2 H 2 O phase at 100 GPa. We also analyzed the electron density around bond critical points (BCPs) ρðr BCP Þ based on the atom-in-molecule (AIM) theory [56]. The density ρðr BCP Þ reflects the relative bond strength, while its Laplacian is used to classify the chemical interactions: Laplacian values are negative for covalent bonds and positive for hydrogen bonds and vdW interactions. These data are available in Table III in the Supplemental Material [48]. We found that hydrogen bond interactions in the Pnma HeNH 3 phase are quite significant while they are small in the I4 1 md He 2 H 2 O phase and almost absent in Fd3m He 2 H 2 O. Therefore, both RDG and BCP analysis show denser and more discrete interactions in Pnma HeNH 3 , which may hinder the diffusion of helium atoms. These indications inspired us to look into crystal structures and diffusion barrier in the following.
It is well known that vdW interactions are weaker than hydrogen bond interactions; therefore, the atoms with pure vdW interactions should diffuse more readily than the hydrogen bonded ones. However, for the superionic state of Pnma HeNH 3 , helium atoms remain almost fixed and the hydrogen atoms become freely diffusive. We therefore analyzed the structural environments of helium atoms in these three crystal structures, as shown in Figs. 5(d)-5(f). In the Pnma HeNH 3 phase, each helium atom is enclosed within a distorted prismatic sublattice connected by N─H covalent bonds and H─N hydrogen bonds, as it roughly occupies the octahedral site of a perovskitelike lattice. Such a compact sublattice can place restrictions on the motions of the helium atoms and can reduce their tendency to diffusion. In the Fd3m He 2 H 2 O compounds at 100 GPa, the water ice sublattice forms hexagonal open channels with side length of 2.22 Å, and helium atoms can diffuse quite easily along these channels. Meanwhile, helium atoms stay relatively far from the ice sublattice, the closest He-H distance is about 1.81 Å. In the Pnma HeNH 3 at 100 GPa the closest He-H distance is about 1.65 Å, which is much shorter. These results are consistent with the conclusions from the charge density analysis. We also constructed the probability density function of the hydrogen atoms based on the atomic trajectories obtained from the simulations at 120 GPa and 2000 K, as shown in Fig. S16 [48], which illustrates the dynamical distribution of hydrogen atoms in the superionic HeNH 3 Pnma phase. The helium atoms are retained in narrow and curved channels along the h100i direction and there is also a high energy barrier between helium atoms, which further constrains the diffusion of the helium atoms. It should be noted that superionic hydrogen will not naturally open diffusion channels for helium since the hydrogen diffusion coefficient is much higher than that of helium: any instantaneous opening in the hydrogen-bonded ammonia network is too short-lived to allow large-scale motion of the helium atoms.

IV. CONCLUSION
We systematically explored the helium ammonia system up to 500 GPa and found three stoichiometries of stable compounds with eight new crystalline phases. Among them the Pnma phase HeNH 3 and P2 1 =c phases of He 2 NH 3 feature a highly distorted perovskitelike structure in which ammonia is arranged in the octahedra, while half of the voids are empty and the other half are doubly occupied in the Pnma phase HeNH 3 and all are doubly occupied in the P2 1 =c phase of He 2 NH 3 . All compounds in the HeðNH 3 Þ 2 stoichiometry belong to a host-guest prototype in which all of the nitrogen atoms arrange in four-and eight-membered rings, and the helium atoms are distributed at the center of the eight-membered rings. In particular, the I4cm phase becomes stable at pressures as low as 3 GPa and its band gap has an abnormal pressure dependence. Furthermore, we carried out AIMD simulations to explore the dynamical properties of HeNH 3 and proposed a new pressure-temperature phase diagram of helium ammonia compounds. This phase diagram features various thermally excited states-plasticity, superionicity, and the full melting line. While ground state phase transitions can be hindered by reaction kinetics [34], thermal excitations are in principle accessible with present-day high-pressure techniques applied to molecular systems up to the megabar regime: rotational disorder in plastic phases through quasielastic neutron scattering [57], superionicity via spectroscopy [27] or conductivity measurements [29], and melting via x-ray diffraction [58].
We found a larger temperature range for superionic hydrogen in helium ammonia compounds than in pure ammonia, shifting the melting line higher and there- fore closer to icy planet isentropes, which may indicate a higher possibility for the existence of helium ammonia compounds under planetary conditions. Whereas near the boundary of the superionic and fluid phases in our helium ammonia phase diagram, the helium diffusion coefficient is much lower than that of hydrogen. Which is quite different from the previously studied helium-water system [13]. This is counterintuitive since the weak interaction between helium and the ammonia network suggests that helium could become mobile more easily than the tighter bound protons. We attribute this phenomenon to the compact ammonia sublattice held by covalent bonds and hydrogen bonds which enclose the helium atoms and form a highenergy barrier that reduces their tendency to diffusion, even if the hydrogen bond network itself rapidly changes as a consequence of the diffusive protons.

A. Crystal structure search
We used a variable-composition structure-prediction algorithm, as implemented in a machine-learning accelerated crystal structure search [59], in which we combined the Bayesian optimization and a Gaussian process model to improve the search efficiency and diversity. This method has been successfully applied in many systems, including helium water compounds [13], T-graphene [60], and metal pentazolate salts [61], etc. We have performed extensive variable-composition searches on different pressures, such as 50, 100, 300, and 500 GPa, as well as fixed stoichiometry (He∶NH 3 ¼ 1∶2, 1∶1, and 2∶1, etc.) at these pressures (not every combination). The crystal structures obtained were cross-checked with results from the ab initio random structure searching approach [62,63], which came up with very similar results.

B. Density functional theory calculations
In this study DFT calculations were performed using the Vienna ab initio simulation package (VASP) [64], accompanied by the projector augmented-wave (PAW) method [65], with the normal version potentials, in which the 1s 2 , 2s 2 2p 3 , and 1s 1 electrons are treated as valence electrons for He, N, and H, and the generalized gradient approximation in the Perdew-Burke-Ernzerhof (PBE) exchange correlation functional [66]. To take the vdW interactions into account, we used DFT-D3 [49] corrections in all static structural calculations. Other charge-density-dependent vdW-DF corrections of Langreth and co-workers [67][68][69] (such as optB88 [70], optPBE [70], and rev-vdW-DF2 [71]), more advanced functionals (SCAN [72], HSE06 [73], PBE0 [50], and RPA [74]), and the hard version of pseudopotentials with smaller core radii are used to crosscheck the results; more details about the methods can be found in the Supplemental Material [48]. All of the structures presented were optimized by employing a plane wave cutoff of 720 eV and dense Monkhorst-Pack Brillouin zone sampling grids with resolutions of 2π × 0.025 Å −1 , leading to ionic and cell optimizations with energy and force convergences better than 10 −6 eV and 0.002 eV=Å, respectively. The CASTEP code [75] was used for the AIRSS searches, and similar convergence parameters as used in VASP were employed. Phonon calculations were performed with the phonopy package [76].
C. Ab initio molecular dynamics simulations AIMD simulations based on the Born-Oppenheimer approximation implemented in VASP [59] were performed in 2 × 3 × 2 supercells for Pnma HeNH 3 (with 240 atoms) and 2 × 2 × 2 supercell for I4cm HeðNH 3 Þ 2 (with 288 atoms), both with Γ-centered k point in the Perdew-Burke-Ernzerhof exchange correlation functional. A normal precision and a cutoff energy of 720 eV were adopted to ensure energy convergence of better than 10 −5 eV. The convergence tests of the k mesh and the size effects are available in Table II in the Supplemental Material [48]. We adopted the canonical NVT ensemble using a Nose-Hoover thermostat [77] with SMASS ¼ 2, lasting for 7 ps with a time step of 1 fs, and we allowed 2 ps for thermalization and then extracted data from the last 5 ps. Some trajectories were extended to 12 ps to check the stability of the simulations.
To obtain reliable diffusion coefficients, we carried out a 22 ps simulations under the same condition, compared with the superionic one under 100 GPa and 2000 K, as shown in Fig. S17 [48]. Also, the vdW corrections do not affect the diffusion behaviors, as shown in Fig. S18 [48].
We calculated the diffusion efficient from the averaged mean-square displacement and cross-checked by velocity autocorrelation function (VACF) by the formula below.
Furthermore, we can extract the vibrational density of states (VDOS) from VACF by Fourier transform: We formally normalized VDOS to the 3 frequency spectrum per species: Z ∞ 0 dωF α ðωÞ ¼ 3: Finally, the nuclear quantum corrections [23] can be expressed by where m α is the mass of atoms for species α and k B is Boltzmann's constant and the total mass of atoms is defined as m ¼ P α N α m α . We specified the time step dτ of 1 fs. We cross-checked the nuclear quantum effect by path integral molecular dynamic with a colored noise method, implemented in the i-PI package [78].

D. Charge density analysis
The atoms-in-molecules (AIM) [56] and reduced density gradient (RDG) analyses [55] of the electron density ρðrÞ were performed by the CRITIC2 code [79].