Hydrogen segregation and its roles in structural stability and metallization: silane under pressure

We present results from first-principles calculations on silane (SiH4) under pressure. We find that a three dimensional P-3 structure becomes the most stable phase above 241 GPa. A prominent structural feature, which separates the P-3 structure from previously observed/predicted SiH4 structures, is that a fraction of hydrogen leaves the Si-H bonding environment and forms segregated H2 units. The H2 units are sparsely populated in the system and intercalated with a polymeric Si-H framework. Calculations of enthalpy of formation suggest that the P-3 structure is against the decomposition into Si-H binaries and/or the elemental crystals. Structural stability of the P-3 structure is attributed to the electron-deficient multicenter Si-H-Si interactions when neighboring silicon atoms are linked together through a common hydrogen atom. Within the multicenter bonds, electrons are delocalized and this leads to a metallic state, possibly also a superconducting state, for SiH4. An interesting outcome of the present study is that the enthalpy sum of SiH4 (P-3 structure) and Si (fcc structure) appears to be lower than the enthalpy of disilane (Si2H6) between 200 and 300 GPa (for all previously predicted crystalline forms of Si2H6), which calls for a revisit of the stability of Si2H6 under high pressure.

T c ~ 16 K to 260 K. Experimental determinations of high-pressure structures of SiH 4 have been carried out using various techniques. Degtyarevaet al. observed from powder x-ray diffraction patterns that SiH 4 remains four-fold to around 25 GPa, and likely transforms to a new phase at higher pressures 9 . Chen et al. observed that SiH 4 sample turns opaque above 27 GPa, which, in combination with measured Raman and IR reflectivity, led to the suggestion that SiH 4 metallizes near 60 GPa 10 . Eremets et al. later reported the metallization of SiH 4 to occur between 50 and 65 GPa based on electrical resistance measurement, and also found this phase to be superconducting(maximum T c ~ 17.5 K) 11 . The crystal structure of the metallic phase (P6 3 ), however, was unusual: its volume shrinks nearly 60% from its precursor the P2 1 /c structure at the phase transition. According to Degtyareva et al., this implies that the observed metallization actually arisesfrom contaminations in the system rather than from the SiH 4 sample 12 . This was supported by the x-ray diffractionstudy of Hanfland et al. 13 in which SiH4 was found to undergo pressure-inducedamorphization at pressures above 60 GPa recrystallizing at 90 GPainto an insulatingI4 1 /astructure original proposed from ab initio calculations by Pickard and Needs 14,15 . At about the same time, Strobel et al. reported a comprehensive experimental study which confirmed Hanfland et al.'s observations 16 . Strobel et al. observed that SiH 4 darkens above ~50 GPa, and attributed it to a partial loss of crystallinity. Above 100 GPa, SiH 4 recrystallizes intothe I4 1 /a (or I-42d) structure, and remainsin this structure to at least 150 GPa.
In addition to pure silane, the mixtures of silane and molecular hydrogen have also been investigated. In mixtures, hydrogen molecules are perturbed by SiH 4 . Such intermolecular interactions can destabilize the covalent bonds in H 2 , showing promise of obtaining a metallic state at moderate pressures. Almost concurrently, Strobel et al. 17 and Wang et al. 18 reported the synthesis of novel SiH 4 -H 2 complexes under high-pressure conditions. Using power X-ray diffraction, Strobel et al. observedthe formation of a stoichiometric SiH 4 (H 2 ) 2 crystal near 7 GPaand identified its structure as a cubic F-43m structure. Using Raman spectroscopy, Wang et al. studied the behavior of SiH 4 -H 2 fluids and established a binary eutectic point at 72 mol% H 2 and 6.1 GPa. Above the eutectic point the fluid crystallizes in a solid solution. The SiH 4 -H 2 phase diagrams revealed in both studies are very similar. Detailed electronic structures and chemical bonding of the SiH 4 -H 2 mixtures have been subsequently examined through the use of first-principles methods [19][20][21][22][23] .
While recent studies revealed a great deal of information on the behavior of SiH 4 under pressure, the jury is still out on the initial proposal. As the research progresses, the reported metallization of SiH 4 near 55 GPa is now facing challenges. Theory suggests that SiH 4 reaches a metallic state at a much greater pressure (220 GPa), through a phase transition to a Pbcn structure 8 . Motivated by the significant interests in this subject, in the present study we employ theoretical structural searches and property predictions to investigate SiH 4 in the post-I4 1 /a region. We systematically searched for stable structures of SiH 4 at pressures above 200 GPa, using a heuristic algorithm based on particle swarm optimization (calypso method) 24,25 . We predict a new structure of SiH 4 above 241 GPa. This structure (P-3) is computed to be metallic and perhaps also superconducting as inferred from the electron-phonon coupling. Our prediction helps to address the phase diagram of SiH 4 at high pressure; while its significance goes beyond. Disilane (Si 2 H 6 ) is the next member in silicon hydride series (Si n H 2n+2 ), which was suggested in theory to become stable under pressure, stable with respect to decompositions into other binaries and/or the elemental crystals 26 . Using similar argument as for SiH 4 , solid Si 2 H 6 was predicted to be metallic and superconducting. With the addition of the P-3 structure in the SiH 4 phase diagram, however, according to our convex hull calculations all previously predicted Si 2 H 6 structures appear to be thermodynamically unstable with respect to the decomposition of Si and SiH 4 to at least 300 GPa. This therefore suggests a revisit of previous studies on solid Si 2 H 6 and a reconsidering of its very stability at high pressure.

Methods
Our structure prediction is based on a global minimum search of the enthalpy surfaces obtained by ab initio calculations at a constant pressure, through CALYPSO (crystal structure analysis by particle swarm optimization) methodology and its same-name code 24,25 . For silane, disilane, silylene and lilicane, structure predictions with upto 6 formula units (f.u.) in the simulation cell were performed at 200 and 300 GPa, respectively, however, for the complex Si 5 H 18 , simulation cells containing only 1 f.u. were considered. The ab initio calculations were performed using density functional theory 27 within the Perdew-Burke-Ernzerhof (PBE) parameterization 28 of the generalized gradient approximation (GGA) as implemented in the Vienna ab initio simulation package (VASP) code 29 . The all-electron projector-augmented wave (PAW) method 30 was adopted with 3s 2 3p 2 and 1s 1 as valence electrons for Si and H atoms, respectively. Planewave energy cutoffs of 500 eV and 700 eV, and uniform Monkhorst-Pack (MP) meshes 31 for Brillouin zone (BZ) sampling with resolutions of 2π × 0.06 Å −1 and 2π × 0.03 Å −1 were employed in the structure predictions and subsequent calculations (e.g. of thermodynamic stability), respectively. The errors in the enthalpy differences of the structures studied in the structure predictions and the subsequent calculations were found to be in the order of 1 meV/atom and less than 1 meV/atom, respectively. The contribution of the electronic entropy to the enthalpy differences aroused by the MP method is less than 0.1 meV/atom, and negligible. The phonon density of states was calculated in harmonic approximation by the finite displacement method 32 as implemented in the PHONOPY program 33 , and the vibrational free energy was estimated therefrom 34

Results and Discussion
In Fig. 1, the calculated enthalpies (H = E + pV) of the candidate structures of SiH 4 are compared between 200 to 300 GPa. In this pressure range, our structural searches successfully predicted a new structure with the P-3 space group, as well as revealed the previously observed/suggested structures, i.e., structures with P2 1 /c, Fdd2, I4 1 /a, C2/c and Pbcnspace groups 8,9,14 . The experimentally identified I4 1 /a structure 11,14,16 is calculated to be the most stable phase at the low-pressure end. The I4 1 /a to Pbcn structural transformation is calculated to occur at about 225 GPa, in a good agreement with previous study 8 . Near 241 GPa, the newly predicted P-3 structure becomes the most stable phase of SiH 4 . Enthalpy of the P-3 structure is constantly lower than that of the C2/c structure 14 , by more than 0.08 eV per formula unit. As well, we predict a metastable C2/m structure that also appears to be more stable than the C2/c structure. The optimized structural parameters of the P-3 and C2/m structures are listed in Table 1. To account for the temperature effects, we estimated the vibrational free energies (at 300 K) for the I4 1 /a, Pbcn and P-3 structures at four different pressures, using a harmonic method 34 (Table 2). As shown  in the inset of Fig. 1, the inclusion of the vibrational free energies does not alter the phase transition sequence, but shift down the two transition points in pressure (see explanation later). It is worth noting that the small proton mass may induce significant quantum nuclear effects, which can deviate the system from a harmonic description, and therefore, shift the two transition pressures further. A rigorous study of this should turn to explicit calculation of the free energy with anharmonic contributions from both the thermal and the quantum nuclear effects includedx 37 . This could be achieved by performing ab initio path-integral molecular dynamics simulations 38 , however, which is beyond the computational load we can afford nowadays. Here, we just speculate that the temperature effects will further stabilize the P-3 structure with respect to the other two structures. Moreover, the appearance of a molecular bonding environment in the P-3 structure (as will be discussed below) suggests that the van der Waals (vdW) interaction may play an important role in the phase transitions. To evaluate this effect, the enthalpies of I4 1 /a, Pbcn and P-3 structureshave been recalculated usingoptB88-vdW functional [39][40][41][42] . The results indicate that the vdW interaction does not change the phase transition sequence, but increase the transition pressures of I4 1 /a → Pbcn and Pbcn → P-3 from 225 and 241GPa to 242 and 285 GPa, respectively. The P-3 structure is shown in Fig. 2(a,b) in two different views. Its unit cell contains 5 SiH 4 formula units. The extended structure can be described as a polymeric Si-H framework intercalated by quasi-hexagonal layers of H 2 units. Isolated H 2 units in Group-IV hydrides are not entirely new -previous studies [43][44][45][46][47][48] suggested them to exist in solid germane (GeH 4 ), stannane (SnH 4 ), and plumbane (PbH 4 ). In SiH 4 , the H 2 units are sparsely populated with negligible intermolecular interactions. This is in contrast to heavier Group-IV hydrides in which the H 2 units are strongly correlated with a tendency for H 2 -H 2 pairing. To analyze the bonding environment, we employed the electron localization function (ELF) 49 as a measure of relative electron localization with respect to a uniform electron gas. The ELF value represents a probability (0 to 1), with large values identifying places in the structure where electrons of opposite spins present, e.g., in cores, bonds, or lone pairs. As Fig. 2(c) shows, the intramolecular regions in the H 2 units have a very large ELF value (close to 1) resulted from a high degree of electron paring. In fact, the ELF map of the H 2 units is consistent to that of a gas-phase H 2 molecule, suggesting that these H 2 units are inact. This conclusion is resonating with the exceptionally short H-H distance in the H 2 units, i.e., 0.75 Å (calculated at 300 GPa), not so different from that of a H 2 molecule.
In the Si-H framework, each silicon atom is coordinated to 12 nearest neighbor hydrogen atoms, which are then bridged to other silicon atoms. Overall, a third of the bridging hydrogen atoms are shared by 4 silicon atoms and the rest are shared by 3, yielding a Si 5 H 18 stoichiometry. Such Si-H geometry goes far beyond the ubiquitous four-fold coordination of Group-IV elements, and on the face of it, one sees the enigma of silicon atoms forming more bonds (12) than the number that seems to be sufficient for their available valence electrons (4). This 'forbidden' increment of crystal coordination, however, is a commonplace in materials under high pressure, and can be explained by the schemes of electron-deficient multicenter bonding 50 . A quantum Atoms in Molecules (AIM) analysis 51 (Table 3) reveals that the bridging hydrogen receive ~0.7 e − per atom from neighboring silicon (which is also seen in Fig. 2(c)). In contrast, the H 2 units almost have no charge transfer from/to the Si-H framework. At this point, silicon atoms have much reduced s/p mixing in its delocalized multicenter bonding, compared with that in molecular SiH 4 (which favors definite coordination geometry, as we know). A similar electron-deficient scenario is also encountered in the other two candidate structures, I4 1 /a and Pbcn, but the difference is that in these two structures the bridging hydrogen atoms are shared by only 2 silicon atoms 8,14 . Each silicon atom is coordinated to 8 hydrogen atoms; pairs of which bridge to four neighboring silicon atoms, forming two three-center-two-electron (3c-2e) bonds. This bonding arrangement is a prototypical electron-deficient geometry that is commonly known for the diborane molecule (B 2 H 6 ) 52 .
Two interesting observations arise from a comparison of the P-3 structure to the I4 1 /a and Pbcnstructures. At first, the P-3 structure has higher coordination geometry which corresponds to a smaller volume. For example, at 250 GPa, the calculated volume for the I4 1 /a, Pbcnand P-3 structuresgradually decreases in a sequence of 13.59 Å 3 /f.u., 13.43 Å 3 /f.u. and 13.17 Å 3 /f.u., respectively. Evidently, the P-3 structure would have the smallest pV work, which provides energetic advantages. Secondly, in the P-3 structure the bridging hydrogen atoms are constrained by more silicon atoms. This yields reduced frequencies for the Si…H…Si bending modes in the mid-frequency range  THz, at 150 GPa), which is the main reason why the vibrational free energy of the P-3 structure is constantly lower than the other two structures (  Table 2. Calculated vibrational free energies (eV/f. u.) of the I4 1 /a, Pbcn, and P-3 structures. △1 (△2) is the energy difference between the P-3 structure and the I4 1 /a (Pbcn) structure.
is mainly resulted from Si lattice modes, so it is very similar for all three structures. The P-3 structure also has high-frequency H-H vibron modes above 100 THz. These modes however do not have notable contributions to the vibrational free energy due to their low density of states. At 150 GPa, for example, the H-H vibrons only contribute 0.04 eV/f.u. to the vibrational free energy of the P-3 structure. An important issue is worthy of investigation at this point. The appearance of the H 2 units in the P-3 structure, as well as the short H-H contact, seems to indicate a trend of segregation toward the Si 5 H 18 + H 2 limit. As crystalline SiH 4 is known to decompose between 50 ~ 100 GPa 16 , it is not unreasonable to speculate it to decompose again at higher pressures. To this end, we examined the phase stability of the P-3 structure with respect to the decompositions into other possible Si-H binaries and H 2 . Four stoichiometries were used for the binaries, namely, Si 5 H 18 , Si 2 H 6 (disilane), SiH 2 (silylene), and SiH (silicane). To find the most stable structure of these stoichiometries, we performed additional structure   Table 3. Bader charges of the P-3 structure calculated at 300 GPa. Negative (positive) sign indicates an electron loss (gain) for the particular atom.
searches at 200 and 300 GPa. The lowest-enthalpy structure that arises from the search was employed as the candidate for each stoichiometry. In Fig. 3, the calculated enthalpy of formation (△ H f , with respect to elemental crystal of H 2 and Si) is shown for all four stoichiometries. Here a negative △ H f means that the binary phase is more stable than the elemental crystals, while the convex hull of the △ H f values (solid lines) connects the stable phases; in this case they are H 2 , Si, and SiH 4 . Clearly, SiH 4 (P-3 structure) is the most stable Si-H stoichiometry at 200 and 300 GPa. The △ H f of the P-3 structure relative to Si 5 H 18 and H 2 (Cmca12 53 ) structure is negative which confirms its stabilities. It is perhaps not surprising that Si 5 H 18 , SiH 2 , and SiH crystals are predicted to be less stable than SiH 4 , since these species are either reactive intermediates or not known to exist experimentally. A significant finding is that the Si 2 H 6 crystal appears to be unstable as well, with respect to the decomposition of Si + SiH 4 crystal. To our best knowledge, Si 2 H 6 is the only known higher silanes that can readily be prepared in laboratory, usually by the reaction of silicon chloride (Si 2 Cl 6 ) with lithium aluminum hydride (LiAlH 4 ). Due to the weaker Si-Si bond (226 kJ mol −1 ), Si 2 H 6 decomposes slowly even at room temperature 54 . On the other hand, previous theoretical studiessuggest that the stability of Si 2 H 6 can be enhanced by applying pressure 26,55 . A series of crystalline Si 2 H 6 polymorphs, with the space group Cmcm, C2/c, and Pm-3m, have been predicted to be thermodynamically stable and superconducting at pressures above 190 GPa. On the contrary, new calculations using the P-3 structure as the most stable structure of SiH 4 show that these Si 2 H 6 polymorphs would all have positive △ H f with respect to the crystal of SiH 4 + Si between 200 and 300 GPa (Table 4). This outcome therefore suggests a revisit on the stability of crystal Si 2 H 6 at high pressures. It should be noted here that the choices of the Si-H binaries as the decomposition products are based upon the known Si-H stoichiometries observed at ambient pressure. It is possible that other Si-H binaries, in additional to what have been considered here, can be stabilized at high pressure.
It is worth noting that another recent theoretical study 56 reports an interesting result on the disproportionation of GeH 4 at high pressure. In this study, zero-point energies were estimated using the   harmonic approximation, which is very similar to the approach adopted in the present study. Under this consideration, a new stoichiometry for the Ge-H binaries, namely GeH 3 , is predicted to become energetically stable with respect to the decomposition of elemental crystals of Ge and H 2 near 175 GPa. The GeH 4 , on the hand, was suggested to be unstable with respect to the decomposition of GeH 3 and H 2 . These results are distinctly different for that of the present study in which the SiH 4 is suggested to be the most stable stoichiometry for the Si-H binaries at high pressures. Experimentally, SiH 4 is the most stable binary in its hydride series (Si n H 2n+2 ) at ambient pressure 54 . At high pressure, SiH 4 decomposes near 50 GPa but recrystallizes above 100 GPa, and remains stable to at least 192 GPa, the highest pressure attempted in the experiments 11,16 . The P-3 structure is predicted to become stable near 177 GPa (with harmonic zero-point corrections), which is well within the experimental stability range for SiH 4 . The study of GeH 3 on the other hand indicates that there are maybe new stoichiometries for this group, perhaps unprecedented at ambient pressures, may become stable at high pressures and this awaits to be discovered in future experiments. Figure 4(a,b) show the calculated electronic band structure and projected density of states (DOS) for the P-3 structure at 300 GPa. It is of considerable interest that this structure is predicted to be metallic. The DOS reveals a pseudogap developed below the Fermi energy (E F ). Similar pseudogaps were also found in calculations on high-pressure GeH 4 and SnH 4 but not on PbH 4 [44][45][46][47][48] . This suggests that the electronic structure of the P-3 structure is not free-electron-like. Since in this case the Fermi level lies in the pseudogap, the DOS at the E F in the P-3 structure is substantially lower than the value in a free electron gas ( / n E 3 2 F , where n is the density of electrons), which limits the possibility of achieving high superconductivity. It perhaps merits mention that the DOS characteristics of the P-3 structure is similar to that of the AlH 3 at high pressure 57 . Since Al has one less electron than Si, the Fermi level of AlH 3 lies in the lower reach of the pseudogap which also results in a low DOS value and ultimately non-superconducting behaviors 58 . In Fig. 4b, hydrogen atoms (mostly bridging ones) contribute to the DOS throughout the entire energy range. This indicates that the hydrogen and silicon electronic states are strongly mixed in the Si-H framework. If the states near the E F can be effectively coupled with the phonon modes, especially the high-frequency ones from hydrogen vibrons, a credible superconductivity can be expected 3 . However, the H 2 units in the P-3 structure are populated sparsely in the lattice so their DOS and phonon DOS are both very low (Figs 4b,d), compared with their counterparts in SnH 4 , GeH 4 and PbH 4 . We therefore do not anticipate exceptional electron-phonon coupling in the P-3 structure at this point. Figure 4(c,d) show the calculated phonon dispersion relations and phonon DOS for the P-3 structure at 300 GPa. The absence of imaginary frequency modes suggests that this structure is mechanically stable. An outstanding feature of the 'H 2 ' intercalation in the structure is that the vibrational modes are divided into subsets in frequency. The H 2 pairs occupy primarily two subsets; one induced by roton modes around 25 THz, and the other induced vibron modes around 118 THz. In these two subsets, the contribution from the Si-H framework is minor, consistent with the fact that the H 2 pairs are inactive. The Si lattice modes dominate the low-frequency subset below 20 THz, while the bridging hydrogen atoms are in the intermediate subset between 30 and 70 THz.
A brief discussion on the possibility of phonon-mediated superconductivity in the P-3 structure is of interest. The methodology was based on an extension of the BCS model in which the attractive, electron-phonon interaction for each phonon mode is treated explicitly. The strength of this interaction is characterized by the electron-phonon coupling parameterλ incorporating the contributions from all participating phonon modes. In Fig. 5, we present the Eliashberg phonon spectral function α 2 F(ω ) and the integrated λ (ω ) ( ) as a function of frequency ω at 300 GPa. The overall integrated λ is 0.63, which is a moderate value for hydrogen-rich materials and comparable to the zero-pressure value for MgB 2 [59][60][61] . The Si lattice modes contribute about 0.23 to the total λ. While the low frequencies of the lattice modes were considered as a disadvantage of attaining superconductivity, their reasonably strong EPC compensates it. The intermediate-frequency H modes contribute 0.36 while the high-frequency vibron modes only contribute 0.04. The logarithmic average of the phonon frequency ω log is obtained as 1320 K. The electron-screened repulsive interaction is represented by the Coulomb pseudopotential μ * . For typical phonon-mediated superconductors, the values of μ * between 0.1-0.13 are generally considered as reasonable 3 . The superconducting critical temperature T c is estimated from the McMillan formula 62 along with the Allen-Dynes correction ( 63 , using an empirical value of μ* = 0.12. The estimated T c of P-3 structure turns out to be 32 K at 300 GPa. A slightly improved estimate of the T c of 35.1 K is made by directly solving the Eliashberg equations 2 with the calculated α 2 F(ω ), using the methodology previously implemented by us 64 . We note that these two estimates of T c are both based on isotropic gap equations and an empirical value of μ*, which may not be adequate depending on the mechanism of the electron-phonon coupling in SiH 4 . A more accurate estimate of the T c may be made using a fully anisotropic treatment of the gap equations, which encourages future study in this direction. Moreover, it should be pointed out that the quantum nuclear effects may deviate the system from a harmonic description and alter the estimated T c . As manifested by a study of AlH 3 , anharmonicity of atomic motions can cause renormalization of the vibrational modes and suppress the superconductivity 58 , On the other hand, however, anharmonic vibrations were found to enhance the electron-phonon matrix elements, as in the case of disordered materials 65 .

Conclusions
We present a theoretical study of high-pressure phase transitions and the metallization of crystalline silane (SiH 4 ). A new polymorph of SiH 4 (space group: P-3) was predicted from structure searches using a heuristic algorithm based on particle swarm optimization (calypso methodology). This new phase becomes thermodynamically stable at pressures above 241 GPa, replacing the previous suggested C2/c structure in the phase diagram. The P-3 structure is calculated to be metallic and also superconducting with an estimated T c of 32 K at 300 GPa. A unique structure feature of the P-3 structure is the presence of H 2 units that are intercalated with the Si-H framework. Electron-deficient multicenter bonding along the Si…H…Si connections results in the delocalization of valence electrons and accounts for the increment of crystal coordination in the Si-H framework. Calculations of enthalpy of formation suggest that the P-3 structure is stable with respective to the decomposition to other Si-H binaries and/or the elemental crystals.