Coexistence of plastic and partially diffusive phases in a helium-methane compound

Abstract Helium and methane are major components of giant icy planets and are abundant in the universe. However, helium is the most inert element in the periodic table and methane is one of the most hydrophobic molecules, thus whether they can react with each other is of fundamental importance. Here, our crystal structure searches and first-principles calculations predict that a He3CH4 compound is stable over a wide range of pressures from 55 to 155 GPa and a HeCH4 compound becomes stable around 105 GPa. As nice examples of pure van der Waals crystals, the insertion of helium atoms changes the original packing of pure methane molecules and also largely hinders the polymerization of methane at higher pressures. After analyzing the diffusive properties during the melting of He3CH4 at high pressure and high temperature, in addition to a plastic methane phase, we have discovered an unusual phase which exhibits coexistence of diffusive helium and plastic methane. In addition, the range of the diffusive behavior within the helium-methane phase diagram is found to be much narrower compared to that of previously predicted helium-water compounds. This may be due to the weaker van der Waals interactions between methane molecules compared to those in helium-water compounds, and that the helium-methane compound melts more easily.


INTRODUCTION
Hydrogen and helium are the most abundant elements in the universe and are significant constituents of the planets in our solar system [1]. Superionicity of hydrogen in ice and ammonia was discovered in early works [2][3][4]. In the superionic state, some atoms form a fixed sublattice while others diffuse as in a liquid. These states have attracted much interest in planetary and high-pressure science for some decades [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19], because ionic mobility affects thermal and electrical conductivity deep inside planets, and therefore their thermal evolution and ability to sustain magnetic fields. For example, it has been reported that the body centered cubic (bcc) phase of superionic ice transforms to a face centered cubic (fcc) phase, which is also superionic [11]; the latter has reportedly been seen in recent shock wave experiments [18,19]. Sun et al. [14] reported a superionic state in close-packed and P2 1 /c phases of ice at higher pressures, and French et al. [15] constructed thermodynamic potentials for the superionic bcc and fcc ice phases and calculated the phase boundary between them.
Superionicity of hydrogen in ammonia and ammonia compounds has also been studied [10,12,13,16,20,21]. Methane is also an important component of giant planets in addition to water and ammonia [17]. Superionic states have not been found in methane, due to the polymerization and release of hydrogen that occurs in methane at high pressures and temperatures [22][23][24][25]. Such polymerization may result in 'diamond rain' in the icy planets [26]. In addition to the superionic states, the plastic phase in molecular crystals such as ammonia, in which protons rotate around the fixed nitrogen atoms, has been reported to emerge at certain pressure and temperature ranges [10]. The plastic phase and rotational motion in water have been studied C The Author(s) 2020. Published by Oxford University Press on behalf of China Science Publishing & Media Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. under high pressure [27][28][29]. Very recently, Li et al. [30] found that colossal barocaloric effects in plastic crystals may have potential applications in solid-state refrigeration technologies. While hydrogen-rich compounds have been studied extensively, this does not hold for heliumcontaining compounds. Traditionally helium is seen as the most inert element because of its closed-shell electronic configuration, but recently helium has been found to react under high pressure with metals [31,32] and ionic compounds [33][34][35]. Liu et al. [33] attributed the reactivity of helium with ionic compounds to the lowering of the Madelung energy between ions arising from the insertion of helium. At high pressures, helium has also been reported to form van der Waals (vdW) compounds with atoms such as neon [36] and covalent molecules such as ammonia [21], water [37][38][39], nitrogen [40][41][42], carbon dioxide [43] and arsenolite [44].

RESEARCH ARTICLE
Superionicity involving helium has rarely been studied. Technically, helium is expected to diffuse as a neutral entity rather than an ionic entity, but a partially diffusive state involving helium has significant implications for thermal conductivities and viscosities. One such example is the helium-iron compound FeHe [32]. The superionic phase of FeHe occurs at pressures above 2 TPa and temperatures higher than 12000 K. Several superionic phases have recently been found in helium-water mixtures, which show novel behaviors such as simultaneous superionicity of hydrogen and helium [39]. However, the possi-bility of compounds in the helium-methane system under high pressures and the nature of their hightemperature behavior are still open questions. This is despite these species forming significant portions of icy planetary atmospheres and mantles, respectively. Their miscibility is then relevant for the atmospheremantle boundary region, which is expected to feature large compositional gradients [45,46]. We have therefore investigated helium-methane mixtures at high pressures and temperatures.

RESULTS
In this work, searches for helium-methane compounds were performed using a machine learning accelerated crystal structure prediction method [47]. Structural optimizations, ab initio molecular dynamics (AIMD) simulations, enthalpies and electronic structures were calculated using the Vienna ab initio simulation (VASP) code [48] with projector augmented-wave potentials and the Perdew-Burke-Ernzerhof exchange-correlation functional (PBE) [49]. Further details of the calculations are provided in the Supplementary Material.
Thermodynamic stabilities of the newly predicted compounds were estimated from their relative enthalpies of formation compared with those of polymorphs of carbon, hydrogen and helium, and C-H compounds predicted in previous works [23,25,[50][51][52]. Using an optB88-vdW functional [53,54] and hard pseudopotentials, we predicted that a new structure with a helium-methane stoichiometry of 3:1 is stable from 55 to 155 GPa (as shown in Fig. 1(a)) which is much wider than the stable pressure region of pure methane molecular crystals. Because methane decomposes above 100 GPa, here we showed a C-H-He phase diagram rather than a He-CH 4 phase diagram. The formation energy of the compound is about 32 meV/f.u. at 110 GPa ( Fig. 1(b)). We have also computed the enthalpy-pressure relations using the PBE and PBE + D3 functionals [55], as shown in the Supplementary Material, all of the results confirmed that He 3 CH 4 is energetically stable. He 3 CH 4 is a molecular crystal with hexagonal space group (SG) P6 3 mc composed of helium and methane molecules. As shown in Fig. 1(c) and (d), helium chains are inserted into the open channels formed by the methane molecules. The packing of methane molecules in He 3 CH 4 is the same as in the experimental hexagonal closed packed (HCP) methane phase at low pressure [56], which is very different from the high-pressure phases of methane with orthorhombic space groups. Another compound, HeCH 4 with SG P2 1 /c, is stable over a narrow range of pressures around 105 GPa. In addition, we have also found several metastable compounds, including HeCH 4 with SG P2 1 , He 2 CH 4 with SG P3 1 m and P2 1 /m, and He 2 (CH 4 ) 3 with SG Cm. These metastable phases are very close to the convex hull. The phase diagrams under different pressures are shown in the Supplementary Material. He 3 CH 4 is an insulator with a large band gap of 8.9 eV at 50 GPa, as shown in Fig. S11 in the Supplementary Material, it is higher than that of pure methane molecular crystals. The phonon band structures of He 3 CH 4 at 0 and 50 GPa shown in the Supplementary Material indicate the dynamical stability of this compound. HeCH 4 is also found to be an insulator with a large band gap and is dynamically stable. The conditions in the upper mantle regions of icy planets reach tens to hundreds of GPa and thousands of K [45]. To obtain the dynamical properties of He 3 CH 4 we performed AIMD simulations within the pressure range 50-200 GPa and the temperature range 500-3000 K. The averaged mean squared displacements (MSD) of hydrogen, helium and carbon atoms were calculated to study phase transitions induced by temperature and pressure. According to the diffusive behaviors of the different atoms, the states of He 3 CH 4 can be divided into four types: the solid phase, a plastic CH 4 phase, a phase of plastic CH 4 plus diffusive He, and the fluid phase.

RESEARCH ARTICLE
We used three representative trajectories to reveal the differences between the states by comparing their MSD and motions, as shown in Fig. 2. All of the simulations start from the relaxed configuration at 150 GPa, which are independent, but with different temperatures. At 1000 K, the atoms are restricted to their equilibrium positions, and exhibit small vibrations. As shown in Fig. 2(a), the MSD of the atoms are extremely small in the simulations, thus He 3 CH 4 remains in the solid phase at this temperature.
When the temperature increases to 1900 K, the compound transforms to the plastic CH 4 phase. As shown in Fig. 2(b), the averaged MSD of H atoms increases rapidly in a short timescale, and stays about 2.0Å 2 afterwards. This is deemed to be a plastic phase when the methane molecules are free to rotate around the carbon atoms with small fluctuations of C-H bond lengths and angles (Fig. 2(e)). Plastic phases have also been reported in pure methane [57], ice [27,58], ammonia [10,12] and helium-ammonia compounds [21]. Meanwhile, the averaged MSD values of He and C atoms are still very small. Using the rigid molecule approximation [57] we have been able to calculate the theoretical MSD of the H atoms. Suppose that the methane molecules are rigid with fixed radius r C H (the length of a C-H bond). For one of the H atoms, its initial position is r 0 (r C H , 0, 0) in spherical coordinates. Since methane molecules rotate freely, the position of the H atom r (r C H , θ, φ) is distributed uniformly over the sphere of radius r C H . Therefore, the analytical MSD is For the plastic CH 4 phase, the resulting MSD is (r − r 0 ) 2 = 2.142Å 2 , which is very close to the value from our AIMD calculations (2.169Å 2 ).
At higher temperatures (2350 K) methane molecules rotate while helium atoms are diffusive (Fig. 2(c) and (f)), which leads to the formation of a superionic-like He state in He 3 Fig. 1(d)) can explain the anisotropy of diffusion coefficients. Heating the superionic phases of He 3 CH 4 eventually leads to melting of the methane sublattice, which gives rise to the fluid phase.
We can further understand the dynamical processes of the He 3 CH 4 structure at different temperatures by using the radial distribution function (RDF). C-H and C-He partial RDFs are shown in Fig. 3(a) and (b) and others are shown in the Supplementary Material. The C-He, He-He and H-He partial RDFs of the partially diffusive phase (2350 K, green lines) are very similar to those of the plastic phase (1900 K, red lines), but are clearly different to those of the fluid phase (2700 K, blue lines). We analyzed some trajectories of the partially diffusive phase and found that the diffusion of helium atoms maintains the order of the helium chains. In some cases of the partially diffusive phase we can see the plateaus in the averaged MSD of He. The performance of the diffusive helium atoms is similar to collective diffusion of ions in lithium battery cathode materials [59,60]. The atoms jump along the chains or hop from one chain to another and therefore the RDFs remain solid-like. However, in the fluid phase, the sublattices of helium and methane have already melted and the RDFs are liquid-like. Finally, the first peak of the C-H RDFs remains essentially unchanged deep into the fluid phase. This indicates that the integrity of the methane molecules is maintained up to the highest temperature.
The superionic state in helium-methane differs from the few other examples of superionic helium [32,39] because plastic CH 4 and diffusive He states coexist. Plastic and superionic phases also appear in the phase diagram of ammonia and ice [10,12,58], but they cannot coexist since both the plastic and diffusive atoms in ammonia and water are protons.
From the results of extensive AIMD simulations, we proposed a phase diagram of He 3 CH 4 between 50 and 200 GPa and below 3000 K (Fig. 4). The colored dots represent independent NVT (Canonical Ensemble) simulations which are classified by their averaged MSDs and RDFs (tests with the NPT (Isothermal-Isobaric Ensemble) confirm these classifications, see the Supplementary Material). The phase boundaries divide the diagram into four regions: solid, plastic CH 4 , plastic CH 4 + diffusive He and fluid phases. The partially diffusive phase appears at pressures above 70 GPa within a narrow range. Compared with helium-water compounds, He 3 CH 4 is dominated by dispersion interactions between methane molecules, which are much weaker than the hydrogen bonds between water molecules. Therefore the ice frameworks in helium-water compounds are more stable than the sublattice of methane molecules.

RESEARCH ARTICLE
There are many open channels in helium-water compounds, which results in a much wider partially diffusive region in helium-water than that in the helium-methane system. To validate the superionicity of helium we have performed AIMD calculations at 150 GPa using hard pseudopotentials. The results are shown in the Supplementary Material and the appearance of the partially diffusive phase at high temperatures is unaffected.
To investigate the interactions in He 3 CH 4 we have applied a real-space analysis to the heliummethane compound. According to the electron localization function (ELF) shown in Fig. 3(d), strong intramolecular covalent bonds persist between carbon and hydrogen atoms while the interactions between methane and helium molecules are of closedshell character. The types of intermolecular interactions in He 3 CH 4 can be determined using the atoms-in-molecules (AIM) theory [61] and reduced density gradient (RDG) analysis [62]. Based on AIM, a topological analysis of the electron density ρ(r) was carried out to compute atomic Bader charges and search for bond critical points (BCPs). We found that the charge transfer in He 3 CH 4 is less than 0.03e between carbon and hydrogen. This is consistent with the fact that the electronegativities of carbon and hydrogen are similar and methane cannot participate in hydrogen bonding. A BCP connecting a pair of atoms provides important information about the bonding between atoms. We show all of the non-equivalent BCPs and their properties in the Supplementary Material. The first and second BCPs represent the bonds inside methane molecules. The electron densities at these BCPs are much larger than others and the negative Laplacian values indicate the concentration of electrons between the carbon and hydrogen atoms. These are typical features of covalent bonding. Other BCPs with small densities and positive Laplacian values are characteristic of dispersion interactions. Furthermore, reduced density gradient analysis has been applied to the electron density. As shown in Fig.  3(c), the spikes at low density, the low-gradient region (inside the red dashed box) reflects the existence of non-covalent interactions [62] in He 3 CH 4 . The distinct spikes are very near zero, indicating that the type is vdW interactions, which are weaker than hydrogen bonds in helium-water [39] and heliumammonia compounds [21].

DISCUSSION
The isotopic effects in hydrogen are important in some cases [63] and we have accounted for them in the helium-methane compound by comparing the vibrational densities of states (VDOS) of He 3 CH 4 and He 3 CD 4 (Supplementary Material). Within the harmonic approximation, isotopic effects lead to changes in the atomic mass and in the phonon frequencies. The frequency of the highest peak in the partial VDOS drops from about 3700 cm −1 for hydrogen to about 2700 cm −1 for deuterium, with a ratio of about 1/ √ 2. Since strong covalent bonds exist between carbon and hydrogen/deuterium, the highest phonon peaks of carbon have the same tendencies as the isotopic effects. In contrast, the frequency region of the helium VDOS does not change substantially, because the vdW interactions between helium and methane are weak. The isotopic effects can also influence dynamical properties, and therefore we conducted AIMD calculations for He 3 CD 4 at 2350 K and compared the VDOS of He 3 CH 4 and He 3 CD 4 at high temperatures. For carbon and hydrogen/deuterium the reduction of the frequencies also appears with the ratio of 1/ √ 2. For helium, the partial VDOS is very similar. This demonstrates that isotopic effects do not substantially affect the superionicity of the heliummethane compound.
In previous studies nuclear quantum effects were considered for ice [15] and mixtures of methane, ammonia and water [64]. The nuclear quantum corrections for He 3 CH 4 at different densities and temperatures are shown in the Supplementary Material. The quantum corrections for the helium-methane compounds are slightly larger than those for ice [15] because of the higher proportion of hydrogen and helium in the compound and the lower temperature region in our simulations.

CONCLUSION
In summary, we have predicted a helium-methane compound (He 3 CH 4 ) that is stable at pressures relevant to upper mantle conditions of icy planets. The inclusion of He atoms highly changes the packing of methane molecules and the stable pressure region of the helium-methane compound is much wider than that of pure methane. He insertion also changes the electronic properties of methane. For example, the band gap of He 3 CH 4 is larger than that of pure methane. Moreover, the phase diagram of He 3 CH 4 has been investigated and a novel phase of coexistence of plastic methane and diffusive helium has been found. The temperature range of the partially diffusive regime in helium-methane is narrower than that in the helium-water system, due to the weaker interactions between the methane molecules, which results in a relatively fragile framework and an easier transition to the fluid state than in the helium-water system. In addition, we observed anisotropy in the He diffusion, which is related to the structure of the compound. We have also analyzed the interactions in He 3 CH 4 and their effects on the phase diagram in comparison with previously predicted helium-water compounds. This work should be helpful in constructing models of icy giant planets, and it would be very instructive to investigate how much the finitetemperature miscibility of helium and methane reflects our ground state results [65].

METHODS
We used a Bayesian-Optimization-based crystal structure prediction method [47] combined with VASP to predict new compounds for the He-CH 4 system. The prediction results have been checked by ab initio random structure searching (AIRSS) [66]. The optB88-vdW functional [53] was employed to account for vdW interactions in the calculations of enthalpies and electronic structures. A basis set energy cutoff of 720 eV was used except for the AIMD calculations, for which a lower cutoff energy of 625 eV was used for reducing the computational cost of extensive AIMD simulations. The Brillouin zone was sampled with a Monkhorst-Pack k-point mesh with a spacing of 2π × 0.03Å −1 . The phonon dispersion curves were calculated using 2 × 2 × 2 supercells with PHONOPY code [67] to validate the dynamical stabilities of the predicted structures. The AIM [61] and RDG [62] analysis of the electron density ρ(r) was performed using the critic2 code [68]. We used orthorhombic supercells containing 256 atoms to perform AIMD simulations for He 3 CH 4 in the NVT ensemble with -centered k-points sampling. The time step of AIMD was set to 1 fs and all simulations were carried out with at least 3000 steps. Some trajectories were extended to more than 10 ps to confirm the stabilities. In addition, to validate our results we employed hard pseudopotentials for hydrogen and carbon and repeated the calculations. The cutoff energy of 1000 eV was used to calculate pressure-composition phase diagrams. A cutoff energy of 910 eV was used for the AIMD simulations. The nuclear quantum corrections of free energies are calculated from the AIMD trajectories using the method proposed in reference [15].