Direct H-He chemical association in superionic FeO2H2He at deep-Earth conditions

abstract Hydrogen and helium are known to play crucial roles in geological and astrophysical environments; however, they are inert toward each other across wide pressure-temperature (P-T) conditions. Given their prominent presence and influence on the formation and evolution of celestial bodies, it is of fundamental interest to explore the nature of interactions between hydrogen and helium. Using an advanced crystal structure search method, we have identified a quaternary compound FeO2H2He stabilized in a wide range of P-T conditions. Ab initio molecular dynamics simulations further reveal a novel superionic state of FeO2H2He hosting liquid-like diffusive hydrogen in the FeO2He sublattice, creating a conducive environment for H-He chemical association, at P-T conditions corresponding to the Earth's lowest mantle regions. To our surprise, this chemically facilitated coalescence of otherwise immiscible molecular species highlights a promising avenue for exploring this long-sought but hitherto unattainable state of matter. This finding raises strong prospects for exotic H-He mixtures inside Earth and possibly also in other astronomical bodies.


INTRODUCTION
The recent discovery of iron peroxide FeO 2 and its reaction with prominent elements and compounds at pressure-temperature (P-T) conditions corresponding to Earth's core-mantle boundary (CMB) [1][2][3][4][5] has opened exciting avenues for exploring new materials [6,7] and properties that have great impact on major geophysical and geochemical phenomena. Of particular significance are FeO 2 -derived minerals that help elucidate distinct pathways for hydrogen and oxygen cycles as well as mechanisms for water transportation deep inside the Earth's interior, which are crucial to understanding key geological events, especially those vital to creating the conditions for life on Earth [8]. Further studies [9,10] revealed mechanisms for structural formation and evolution of a dense phase of FeO 2 H at CMB P-T conditions. Moreover, a recent work [7] predicted a rare helium-bearing compound FeO 2 He that is robust and viable in Earth's CMB region, offering a compelling scenario for elucidating the origin of enigmatic deep-Earth reservoirs for primordial helium inferred by geochemical evidence [11][12][13][14]. The unique ability of FeO 2 to bond with hydrogen and helium and form new compounds has major implications for understanding the evolution and dynamics of Earth and other giant solar and extrasolar planets.
Hydrogen and helium are the most abundant elements in the universe and play crucial roles in geological and astrophysical environments, but they are known to be inert toward each other across wide P-T and concentration ranges and remain largely immiscible up to multi-megabar pressures and 3000-4000 K temperatures [15][16][17][18][19][20][21][22]. Given their prominent presence and influence on the formation and evolution of celestial bodies, it is of great interest and significance to explore and decipher the nature of interactions between hydrogen and helium, especially possible chemical association that would have considerable impacts in many scientific fields, from chemistry, physics, geoscience to astrophysics [23][24][25]. Recent studies have identified compounds that host both hydrogen and helium, such as in highly compressed helium-water compounds [26][27][28][29], but there is no evidence of direct chemical association of these two elements. A recent experimental study [30]  association of the constituent molecules at pressures < 75 GPa; but later experiments presented opposing conclusions [31,32], showing no evidence for chemical association in H 2 -He mixtures up to 250 GPa at 300 K and pointing to N 2 contamination as the source for the previously observed phenomena [30]. The quest to seek possible reactivity of hydrogen and helium remains an intriguing open challenge.
In this work, we present evidence from firstprinciples calculations showing direct and prevalent chemical association of hydrogen and helium facilitated by their reaction with iron peroxide FeO 2 in forming a rare quaternary compound FeO 2 H 2 He. We conducted extensive crystal structure searches and evaluated the Gibbs free energy of the three reactants and the resulting product to determine their relative phase stability; we also performed ab initio molecular dynamics (AIMD) simulations to assess temperature-driven partial melting of the compound. The results show that FeO 2 H 2 He is viable in a large region of the P-T phase diagram. Most interestingly, in a wide swath of the phase space corresponding to Earth's lowestmantle regions, this quaternary compound stays in a superionic state hosting liquid-like hydrogen inside the FeO 2 He sublattice that remains intact in crystal form. This exotic solid-liquid mixture state of matter makes an unusually conducive environment, promoting close coalescence of hydrogen and helium. An analysis of the H-He radial distribution function (RDF) and bonding features in superionic FeO 2 H 2 He reveals significant direct H-He chemical association as indicated by RDF peaks at short distances with considerable bonding as revealed by interatomic charge characters. These results highlight a compelling case of H-He chemical association, which may be harbored in deep-Earth regions, and the constructed phase diagram provides crucial guidance for exploring solid and superionic phases of FeO 2 H 2 He in laboratory experiments and for modeling interiors of giant solar and extrasolar planets. The computational approach adopted in this work in search of stable quaternary compounds represents the state of the art in the very active and diverse research area of crystal structure prediction, and our reported results may help resolve a significant scientific problem of long-standing and broad interest in many fields of chemistry, physics, geophysical and planetary sciences among others.

RESULTS AND DISCUSSIONS
Our first-principles crystal structure search identified a quaternary compound with stable stoichiometry, FeO 2 H 2 He. Calculated enthalpy data  Table 1) and electronic structures (see Supplementary Fig. 1) are provided in the Supplemental data. While these predicted structures have not been checked exhaustively against all possible decompositions because of computational constraints, they have been extensively examined to assess their viability. A full search process has established that this compound is the most stable in the Fe-O-H-He quaternary compositional space and, crucially, energetically favorable compared to the expected reactants, as indicated in Fig. 1(a). Moreover, we performed phonon calculations, and the results showed that there are no imaginary modes in the entire Brillouin zone from 100 GPa to 300 GPa (see Supplementary  Fig. 2), thus confirming the dynamical stability of the predicted FeO 2 H 2 He structures. Results in Fig. 1 show that the sublattice of the Fe-O-H framework in FeO 2 H 2 He is identical to that of the P-3m1 phase of  Fig. 3), but does not change the main conclusion of this work.
In assessing fundamental behaviors of the newly identified FeO 2 H 2 He crystal phases, a crucial task is to establish their P-T phase diagram. To this end, we performed two sets of calculations. First, we calculated phonon dispersions and the corresponding phonon density of states (PDOS) within the harmonic approximation as implemented in the PHONOPY code [37], and then used the obtained PDOS as input to evaluate the vibrational contribution to the entropy of each phase. Combining with the total internal energy, pressure and volume determined by first-principle calculations, we obtained for both the R-3m and Pnnm phases their Gibbs free energy, which allow determination of the phase boundary as shown in Fig. 2. Second, we performed AIMD simulations for the two FeO 2 H 2 He phases in the pressure range of 70-300 GPa over a temperature range of 0-6000 K to evaluate thermodynamic stability and temperature driven diffusive and melting behaviors by examining mean square displacements (MSDs) of the Fe, O, H and He atoms in the crystal structures. Some representative MSD data are presented in Fig. 3 to showcase the typical behaviors and illustrate the assessment process and criteria. Following this procedure that has been applied to sufficient P-T points, we obtained the boundaries, as indicated in the phase diagram shown in Fig. 2, between the three distinct states: (i) solid, where all atoms remain near their equilibrium crystal lattice sites, (ii) diffusive H, where H atoms show significant deviations from their equilibrium sites while all other atoms remain near their equilibrium sites, and (iii) fluid, where all atoms deviate from their equilibrium sites. In all the simulated cases, increasing temperature first drives H atoms off their equilibrium sites into a diffusive phase, followed by diffusive behaviors of all other atoms that happen almost concurrently, signifying the melting of the crystal structure into the fluid phase. It is noted that the predicted FeO 2 H 2 He structures may coexist with the dehydrogenation products FeO 2 H, H 2 and He toward the low-pressure side of the P-T phase diagram shown here; however, results from AIMD simulations show that the FeO 2 H 2 He phases are at least metastable and there is no immediate phase transition to another FeO 2 H 2 He structure in the P-T space studied here as indicated by the stable atomic positions (see below). We adopted a Boltzmann distribution description for phase coexistence probability [38] of FeO 2 H 2 He relative to its dehydrogenation products at selected pressures along the geotherm [39] (see Fig. 2), and the results indicate robust presence of FeO 2 H 2 He in the P-T regions of interest.
We show that each of the solid, diffusive H, and fluid states of the two FeO 2 H 2 He phases has a large stability field in the P-T phase diagram (Fig. 2), determined by extensive energetic (Gibbs free energy), dynamic (phonon dispersion) and thermodynamic (AIMD) calculations. It is noted that the stability field of the diffusive H state of the R-3m phase includes the Earth's geotherm in the vast lowermantle region, covering the core-mantle boundary (CMB), indicating that the superionic state of this extraordinary quaternary compound may exist in Earth's lower-mantle regions. Moreover, this superionic state, in which H atoms move in a liquidlike state in the still intact solid Fe-O-He sublattice, greatly enhances the probability of bonding interactions of H atoms with other atoms in the compound, especially the long-sought H-He chemical association. To explore this intriguing prospect, we performed extensive computational studies and systematic analysis of the H-He bonding in FeO 2 H 2 He under pertinent P-T conditions. We calculated RDFs from AIMD simulation data to assess interatomic distances in the solid and superionic phases, along with an examination of select short bond lengths. We also performed a systematic evaluation of the electronic density and the Laplacian of the electronic density at the bond critical point, which is the saddle point along neighboring bond paths for interatomic charge distribution commonly used in characterization of bonding in the context of QTAIM [40].
We first set a key reference by examining the bonding and charge states of the two solid FeO 2 H 2 He phases at representative pressure points at 0 K. The H-He bond length in the R-3m phase at 100 GPa is 1.67Å and the value for the Pnnm phase at 300 GPa is 1.50Å. While these bonds are considerably shorter than the sum of the van der Waals radii for H and He atom (2.50Å), a QTAIM analysis (see results in Table 1 and discussions below) showed no evidence of H-He bonding. Calculated RDFs (Fig. 4) exhibited sharp peaks reflecting well defined and separated interatomic distances in these solids.
At elevated temperatures and pressures pertinent to deep-Earth regions, the RDFs show considerable thermal fluctuations (Fig. 4), indicating greatly enhanced ranges of atomic movements. While all the RDFs are thermally broadened, those involving H exhibit more pronounced changes because of the diffusive nature of the H atoms in the superionic phases of FeO 2 H 2 He. In particular, the H-He RDFs undergo the most dramatic broadening, with substantial shifts of the H-He RDF weights to the lower distances. Such remarkably reduced interatomic  separations suggest higher probability for H-He bonding interactions. To verify this scenario, we examined the AIMD simulation data and found significant bonding features in the two superionic phases when the H-He bond length is shorter than 1.3Å, as indicated by the electronic charge density ρ and the Laplacian of the charge density ∇ 2 ρ at the bond critical point where key bonding characters are evaluated in QTAIM to probe the formation of interatomic bonding [40]. In Table 1, we present pertinent results that showcase appreciable H-He bonding comparable to conventional hydrogen bonding [41], which offers compelling quantitative evidence for direct H-He chemical association in superionic FeO 2 H 2 He. Our analysis reveals that the H-He pairs with clear bonding features comprise a large fraction (up to 20-40%) of the total number of H-He pairs in the AIMD simulation cells that fall under the first RDF peaks for the respective superionic phases. This phenomenon indicates a strong tendency for H-He bonding formation in FeO 2 H 2 He at P-T conditions pertinent to deep-Earth environments.
It is noted that the charge analysis performed on a snapshot taken from the MD simulation can be considered representative of the system in a time averaged sense. The comparison with the stable hydrogen bonding is meant to provide a relative reference and measure with an established bonding configuration, with the understanding that the charge configurations in FeO 2 H 2 He are dynamic and any given specific 'bonding' is only transient in nature. On this basis, the present analysis of the dynamic bonding provides a useful description of the superionic state of FeO 2 H 2 He, especially the degree of H-He interaction in such a novel state over a wide swath of the P-T space. We performed calculations for insights into the dynamic charge distribution in FeO 2 H 2 He. The calculated Crystal Orbital Hamilton Populations (COHPs) (see Supplementary Fig. 4) show that FeO 2 H 2 He at 147 GPa and 3100 K host anti-bonding states, and there are overlaps between the H 1s and He 1s states, indicating charge transfer from the H to He atoms. Moreover, we analyzed the data from the molecular dynamics simulations on the nearest H-He distance of 0.9-1.3Å, and the results (see Supplementary Fig. 5) indicate that the lifetime for such bonding states is around 0.01 ps.

CONCLUSION
In summary, we identified a novel quaternary compound FeO 2 H 2 He using an advanced crystal structure search approach and established its P-T phase diagram from extensive first-principles calculations of structural, phono, and molecular dynamics. This compound adopts three distinct phases: solid, superionic with diffusive H and liquid. Each of them occupies a large region in the phase diagram. The superionic phase is viable under the conditions in Earth's vast lower mantle, which includes the enigmatic CMB region where pertinent reactants have been predicted to be present. The superionic nature of this new compound promotes unprecedented close coalescence of hydrogen and helium in a chemically facilitated environment. These findings provide crucial knowledge for deciphering the long-sought H-He chemical association in a large swath of P-T phase space previously thought to be off-limits for H-He miscibility, offering new insights for elucidating novel physics and chemistry of the predicted compound FeO 2 H 2 He in a wide range of high P-T conditions and also for modeling solar and extrasolar giant gas planets.

Structural predictions
Crystal-structure searches for quaternary compounds are computationally demanding; our task here is facilitated by knowing the constituent reactants, namely iron peroxide, hydrogen and helium, which have been predicted to exist in Earth's lowest-mantle regions through geophysical and geochemical processes [1][2][3][4][5]7,[11][12][13][14]. Accordingly, we have explored Fe-O-H-He compositional space with energetic evaluations in relation to FeO 2 , H 2 and He to determine relative structural phase stability.
Our structure search is based on global optimization of potential-energy surfaces using the CALYPSO methodology [42,43], which has been successfully employed in predicting a large variety of crystal structures [6,44,45]. Evolutionary variablecell calculations were performed at 100, 200 and 300 GPa with 1, 2, 3, and 4 FeO 2 H 2 He formula units per cell, retaining 60% lowest-enthalpy structures to produce the next-generation structures by a particle swarm optimization procedure and generating the remaining 40% structures randomly within the symmetry constraint. The number of atoms in the unit cell is in the range of 20-30 for most of the reported structure searching simulations. In light of several Fe-O compounds predicted at 100 GPa and 300 GPa [1,46] Table 2 and Supplementary Table 3 for details). Most searches converge in 30-40 generations with about 1000 structures generated; in all, we have examined about 227 000 structures in this study.

Ab initio calculations
First-principles total-energy and electronic property calculations were carried out using density functional theory (DFT) as implemented in the VASP code [47], adopting the frozen-core all-electron projector-augmented wave (PAW) method [48], with 3s 2 3p 6 3d 7 4s 1 , 2s 2 2p 4 , 1s 1 and 1s 2 treated as valence electrons for Fe, O, H and He, respectively, and the Perdew-Burke-Ernzerhof (PBE) exchangecorrelation functional in the generalized gradient approximation (GGA) [49,50]. Correlation effects among the Fe 3d electrons were treated in the GGA + U approach [51,52], adopting the recently proposed on-site Coulomb interaction U = 5.0 eV and a Hund's coupling J = 0.8 eV [53]. A cutoff energy of 1000 eV for the plane-wave expansion and fine Monkhorst-Pack k-meshes [54] were chosen to ensure enthalpy convergence of several meV/atom.

ZPE, phase-diagram and phonon calculations
To properly describe the quantum effect in hydrogen and helium, we included the zero-point energy (ZPE) in all energetic calculations. To determine the phase boundaries between the solid phases, we performed calculations to determine the Gibbs free energy G = U + PV − TS, where U, P, V, T and S are the internal energy, pressure, volume, temperature and entropy, respectively. The vibrational energy and entropy were obtained from lattice dynamic calculations using the quasi-harmonic approximation as implemented in the PHONOPY code. Phonon calculations were carried out using a supercell approach [55] as implemented in PHONOPY code [37].

Bonding and molecular dynamics calculations
We performed bonding analysis using the quantum theory of atoms in molecules (QTAIM) approach [40]. We also performed ab initio molecular dynamics (AIMD) simulations in the NVT (N-particle, Vvolume, T-temperature) ensemble implemented in the VASP code to probe temperature-driven diffusive and melting behaviors. We employed a simulation time of 10 ps with a time step of 0.5 fs and a 3 times supercell containing 162 atoms in the R-3m phase and a 2 times supercell containing 144 atoms in the Pnnm phase.