Introduction

The unexpected discovery of high-temperature superconductivity in high-pressure hydrides has revamped the hopes of ending the century-long quest for superconductivity at ambient conditions1,2,3,4,5,6,7,8,9,10,11,12,13. In less than five years, superhydrides have established higher and higher records for critical temperatures, starting with SH3 (203 K)3, LaH10 (265 K)4,5, and C-S-H (288 K)9.

While it is very exciting to strive for new superconductors with even higher Tc’s attaining hot superconductivity14,15,16, lowering the required stabilization pressures, while retaining Tc above the boiling point of nitrogen (77 K) is of even greater importance16,17,18,19,20,21,22. In fact, the discovery of a conventional (s-wave) superconductor with these properties would open up a wide array of technological applications in key strategical sectors such as energy conservation, climate change, and medicine. After exhausting the search space of binary hydrides, the focus of superhydride research is rapidly shifting to ternary hydrides, where the parameter space is much larger6,17. We have recently predicted, for example, that a ternary hydride with LaBH8 composition could form an \(Fm\bar{3}m\) ternary sodalite-like clathrate (SLC) structure and remain stable down to a critical pressure pc ~ 35 GPa, with a Tc = 126 K16,19. Our findings were later confirmed by independent studies on the La-B-H system21,23. The predicted pc of LaBH8 is a factor of four lower than in binary hydrides, where pc’s are in the Megabar range, but still too high to envision any large-scale applications for this particular compound. However, through the identification of the \(Fm\bar{3}m\) XYH8 structural template, the discovery of LaBH8 paved the road to the study of a whole new family of potential high-Tc ternary hydrides, where pc and Tc may be improved even further.

In this work, using first-principles methods for crystal structure prediction and superconductivity, we identify two high-Tc alkaline-earth/silicon superhydrides, BaSiH8 and SrSiH8. We predict that both compounds will spontaneously form in the \(Fm\bar{3}m\) XYH8 structure at high pressures (130 and 174 GPa, respectively), and remain dynamically stable down to much lower pressures (pc = 3 and 27 GPa), with superconducting Tc’s of 71 and 126 K, respectively.

We also succeed in deriving a rigorous limit for the stability of BaSiH8, calculating explicitly the energy barrier protecting the metastable \(Fm\bar{3}m\) structure from decomposition (kinetic stability) as a function of pressure, using the Variable-Cell Nudged Elastic Band (VCNEB) method24. We find that, indeed, \(Fm\bar{3}m\) BaSiH8 should remain metastable at pressures well below 130 GPa. However, the kinetic critical pressure pkin determined by the energy barrier is significantly higher (30 GPa) than the value estimated from anharmonic lattice dynamics. Similar discrepancies between dynamical and kinetic pressure of stability may explain the systematic underestimation of the predicted pressures of (meta)stability found in other high-pressure hydrides compared to experiments11,25,26,27.

Results

Initial screening

BaSiH8 and SrSiH8 were first identified through a high-throughput screening of possible substitutions of La and B by neighbouring elements into the \(Fm\bar{3}m\) XYH8 structure, as shown in Fig. 1a. In this structure – Fig. 1b – hydrogen atoms form rhombicuboctahedral cages around lanthanum; boron, being much smaller, fills the six cubic voids surrounding the cages. This structure permits a very efficient realization of a mechanical precompression mechanism observed in SLC binary hydrides25,28. To elaborate: in addition to the large central atom (La), also the smaller atom in the interstitials (B) exerts an additional pressure on the hydrogen sublattice, lowering the minimum stabilization pressure pc19.

Fig. 1: Elemental composition and phase space search.
figure 1

a Dynamical stability of the XYH8 structure at 50 GPa for various combinations of X and Y elements. Dynamically stable and unstable compounds are indicated by red crosses and green ticks, respectively. b Crystal structure of \(Fm\bar{3}m\) XYH8 in the fcc conventional unit cell; H atoms are shown in blue, H-H bonds in grey, and the lighter (heavier) atom in orange (green). c Ternary convex hull for BaSiH8 at 100 GPa - side view; sampled structures are shown as orange dots. d Ternary convex hull for BaSiH8 at 100 GPa - top view. In c and d the color scale (from orange to blue) represents the depth of the hull with respect to the elemental composition, from 0 to −0.7 eV/atom; stable compositions are shown as blue dots and are labeled.

Full size image

We assumed that, considering different combinations of large (X) and small (Y) atoms, the superconducting properties of LaBH8 could be improved even further. Aiming at identifying compounds with low critical pressures, we performed structural relaxations at 50 GPa for all combinations of X = [K, Rb, Cs, Ca, Sr, Ba, Sc, Y, and La] and Y = [B, Al, Ga, C, Si, and Ge] elements in the \(Fm\bar{3}m\) XYH8 template, and evaluated the dynamical stability of the resulting compounds by calculating the harmonic phonon dispersions on a 4 × 4 × 4 grid in reciprocal space.

As shown in Supplementary Fig. 6 of the Supplementary Material (SM), the lattice constants of the relaxed structures exhibit an almost perfect linear dependence on the sum of the X and Y empirical atomic radii R = RX + RY29, with slopes determined by the total number of valence electrons Ne = NX + NY. Out of the 54 compounds investigated, only three compounds with Ne = 6 – green ticks in Fig. 1a – are dynamically stable at 50 GPa: Apart from the already reported LaBH819, we also found two silicon hydrides, BaSiH8 and SrSiH8. Both compounds exhibit a larger R than LaBH8, suggesting that mechanical precompression in both compounds will be more efficient, and hence their pc may be lower. In particular, BaSiH8, where R is considerably larger than in the other two compounds, looks very promising in this respect. Indeed, our harmonic phonon calculations yield pc = 5 and 30 GPa for BaSiH8 and SrSiH8, respectively, lower than the pc = 40 GPa in LaBH8.

Convex hulls and phase diagrams

Having established that BaSiH8 and SrSiH8 are dynamically stable in the \(Fm\bar{3}m\) structure close to ambient pressure, we determined the pressures at which the two compounds may spontaneously crystallize in the \(Fm\bar{3}m\) structure starting from appropriate precursors. To this end, we computed the full ternary convex hulls for the two ternary X-Si-H systems at ambient pressure and at 100 GPa, employing ab-initio variable-composition evolutionary crystal structure prediction methods as implemented in USPEX30,31. To construct each hull, we sampled around 20000 structures and 5000 compositions, including the contribution of the zero-point energy (ZPE) from the ionic vibrations for structures close to the hull. Our results are in good agreement with previously reported structural searches on the Ba-H and Sr-H binary hydrides32,33. Further computational details are given in Supplementary Method 2 and Supplementary Note 2 in the SM. An example is illustrated in Fig. 1c, d, where we show the ternary convex hull of the Ba-Si-H system at 100 GPa. (The analogous convex hull for the Sr-Si-H system at 0 and 100 GPa is shown in Supplementary Figs. 9 and 10 of the SM.)

In the Ba-Si-H system, the BaSiH8 composition is thermodynamically stable, i.e., it lies on the hull at 100 GPa in a distorted P1 phase. The \(Fm\bar{3}m\) phase is slightly higher in enthalpy (32 meV/atom), but becomes enthalpically favourable above 130 GPa. In the Sr-Si-H system at 100 GPa, according to our calculations, SrSiH8 should decompose into SrSiH6 and SrSiH12, which are the stable ternary compositions along with Sr2SiH10. The corresponding structures are shown in Supplementary Fig. 12 of the SM. However, the \(Fm\bar{3}m\) phase lies only 54 meV/atom above the convex hull. From a comparison of the formation enthalpies, we predict that it should become stable at 174 GPa.

Anharmonic lattice dynamics

Preliminary phonon calculations in the harmonic approximation indicated that both BaSiH8 and SrSiH8 experience remarkable phonon softening at lower pressures. More specifically, we find soft phonon modes, mainly at Γ (T2g, Eg, and A1g) and X (Eg, and additionally A2g in SrSiH8), as shown in Fig. 2. Eventually, imaginary frequencies appear for the Eg mode at X below 5 GPa for BaSiH8 and 40 GPa for SrSiH8. (A full set of phonon dispersions at all pressures considered is provided in the SM.)

Fig. 2: Vibrational structure and electron-phonon properties.
figure 2

a Phonon dispersion of BaSiH8 at 5 GPa; the thickness of the dark blue lines is proportional to the mode- and wave-vector resolved coupling constant λq,ν. Solid (dotted) lines correspond to anharmonic (harmonic) results. b Phonon DOS shown as solid black line, α2F(ω) as filled curve, and λ(ω) as solid blue line. Panels c and d as a and b but for SrSiH8 at 30 GPa.

Full size image

Soft-mode behaviour, associated with strong anharmonicity, has been reported for many hydrogen-rich materials25,26,34,35,36. Anharmonic lattice effects have been shown to crucially affect the range of dynamical stability, phonon frequencies and eigenvectors of superhydrides, and their inclusion is essential to obtain accurate estimates of these quantities.

In order to account for anharmonic effects on the phonon dispersions of BaSiH8 and SrSiH8, we evaluated the adiabatic potential energy surface (APES) for every soft mode of a 2 × 2 × 2 wave-vector grid, which includes explicitly the special points Γ, L, and X, and solved the resulting Schrödinger equations. (More information can be found in Supplementary Method 4 of the SM, and Ref. 37).

The harmonic and the anharmonic phonon dispersions for BaSiH8 and SrSiH8 close to their critical pressures of stability are shown as dashed and full lines in Fig. 2. Anharmonicity causes a considerable hardening of the T2g (at Γ) and Eg (at Γ and X) modes, leading to a decrease of the critical pressure of dynamical stability pc: Taking this hardening into account, SrSiH8 is stable down to a pressure of 27 GPa and, even more excitingly, BaSiH8 down to 3 GPa.

Electronic structure and superconductivity

Having determined that both BaSiH8 and SrSiH8 remain dynamically stable close to ambient pressure, we are left with the question whether at these relatively low pressures they can still be considered superhydrides. A first positive indication comes from the analysis of their electronic structure.

In Fig. 3 we show the electronic dispersions and DOS for BaSiH8 (left) and SrSiH8 (right), at 5 and 30 GPa, respectively. Despite a scaling of the total bandwidth due to the different pressures, the band structures are qualitatively very similar. The Fermi level is located just above a large DOS shoulder; in this region and up to −5.5 eV for BaSiH8 (−3.5 meV for SrSiH8), the bands are of purely hydrogen character. This is an important prerequisite for high-Tc conventional superconductivity in hydrides, since it can imply a strong electron-phonon (ep) coupling between electronic states to the light hydrogen sublattice38.

Fig. 3: Electronic band structure and density of states.
figure 3

a Electronic band structure of BaSiH8 at 5 GPa and b corresponding total DOS (solid black) and projected H DOS (filled blue). Panels c and d as a and b but for SrSiH8 at 30 GPa.

Full size image

Similar conclusions can be inferred from the electron localization function (ELF) of the two compounds (see Supplementary Fig. 15 in the SM), where, in line with what we observed in LaBH819, we find an increased charge localization in the vicinity of the H and Si atoms and between H-H, but not between Sr/Ba and H or Si and H. This supports the idea that neither Sr/Ba, nor Si form bonds with H, and thus act on the H sublattice essentially as mechanical spacer, as in binary sodalite-like clathrate hydrides6,25,34,39.

The Fermi surface topology is the same in the two compounds: a large, spheric-like electron pocket is centered around the Brillouin zone center, while a more complicated hole-like network extends around the faces of the Brillouin zone, enclosing the X and W points (see the insets of Fig. 4).

Fig. 4: Superconducting gap as a function of temperature.
figure 4

Δ(T) for BaSiH8 (blue) at 5 GPa and for SrSiH8 (green) at 30 GPa calculated using the fully anisotropic ME theory; dashed lines indicate the mean value as a guide to the eye. The insets in the top right and bottom left corner show the Fermi surface distribution of Δk at 10 K for BaSiH8 and SrSiH8 respectively; blue/red correspond to min(Δ)/max(Δ).

Full size image

These qualitative electronic structure arguments are confirmed by actual calculations of the superconducting properties. Fig. 4 shows the energy distribution of the superconducting gap Δ as a function of temperature T for BaSiH8 at 5 GPa pressure and for SrSiH8 at 30 GPa, obtained by solving the anisotropic Migdal-Eliashberg (ME) equations on a 30 × 30 × 30 k − and q-grid using the anharmonically corrected phonon dispersions with the EPW code40,41. (The superconducting properties at all other considered pressures and further computational details are provided in the SM.) We observe two distinct superconducting gaps: The inset of the figure shows that large Δ values correspond to the Γ-centered, spherical electron pocket, while lower values occur on the hole-like tubular network around the X and W points. The superconducting critical temperatures are predicted to be 71 and 126 K for BaSiH8 and SrSiH8, respectively.

Further details on the origin of the remarkable Tc’s of BaSiH8 and SrSiH8 can be obtained from an analysis of the distribution of their ep coupling over phonon branches. In panel b of Fig. 2, the mode and wave-vector resolved ep coupling λq,ν are overlayed onto the phonon dispersions as fat bands; panel d of the same figure shows the Eliashberg spectral function α2F(ω), and the total frequency-dependent ep coupling parameter λ(ω).

The figure clearly shows that the largest contributions to the total coupling come from low-energy modes; a substantial fraction is associated to soft phonon modes around Γ and X. In the case of BaSiH8, for example, we estimate that the phonons related to Γ − T2g contribute roughly 40% to the total λ and around X − Eg about 35%. These modes are of purely hydrogen character and their patterns can be more easily visualized in terms of the H cubes around Si: The T2g mode at Γ is a centrosymmetric stretching and squeezing along one of the space diagonals of the H cube, while the Eg mode at X is a nondeforming rotation around a face normal of the H cube, with alternating phase in two neighbouring cells.

Figure 5 shows the evolution of the superconducting properties of BaSiH8 and SrSiH8 as a function of pressure from the lowest dynamical stability pressure up to 100 GPa; open circles and solid lines correspond to anharmonic results, while dashed lines refer to harmonic values. Interestingly, Tc is approximately constant with pressure – panel a; this is true even close to pc, where one would expect a strong deviation due to anharmonicity. An inspection of panels b and c shows that this happens because of the compensating effect of phonon hardening on ωlog and λ. The effect is more marked in SrSiH8 and in BaSiH8 close to the instability point, where harmonic and anharmonic dispersions diverge the most. For SrSiH8 at 30 GPa, we find an anharmonic (harmonic) λ of 3.21 (3.96) and an ωlog of 33.1 meV (22.3 meV), and for BaSiH8 at 5 GPa, we find a λ of 1.61 (2.55) and an ωlog of 33.3 meV (23.7 meV).

Fig. 5: Superconducting properties as a function of pressure.
figure 5

a Tc from anisotropic ME theory, b ωlog, and c λ (bottom panel) as functions of pressure for BaSiH8 (blue) and SrSiH8 (green). Circles and solid lines (crosses and dashed lines) indicate anharmonically corrected (harmonic) calculations. The contributions to α2F(ω) from the Eg mode at X are set to zero in the harmonic calculations for SrSiH8 at 30 and 40 GPa, as the corresponding phonon frequencies are imaginary.

Full size image

Metastability

In all theoretical studies of high-pressure hydrides performed so far, the range of metastability of high-pressure phases has been assumed to coincide with the range of dynamical stability, eventually including quantum corrections to lattice dynamics. However, a comparison of these predictions and available experimental data suggests that theoretical estimates of the critical pressure are systematically lower than experimental values11,25,26,27.

Indeed, dynamical stability is only a prerequisite for thermodynamic metastability. The latter is determined by the existence of an energy (enthalpy) barrier protecting a metastable phase from decomposition into other phases42. Attempts to quantitatively estimate the barrier height are extremely rare. For \(Fm\bar{3}m\) BaSiH8, where the dynamical instability pressure pc is extremely close to ambient pressure, this issue is obviously crucial, since the presence of a sufficiently high barrier could be considered the definitive proof of synthesizability.

In this work, we estimated the barrier height using the VCNEB as implemented in the USPEX code24. In this method, a number of intermediate structures (images) are created between the metastable structure of interest and the ground-state structure at a given pressure; elastic forces between each image are added to the ‘physical’ forces, and the whole chain of images thus created is finally relaxed to obtain the energy/enthalpy profile of the transition.

We ran VCNEB simulations for BaSiH8 at six different pressures: 10, 25, 35, 50, 100, and 200 GPa. As end-members for the VCNEB path, we chose the \(Fm\bar{3}m\) phase and a P1 phase, identified as the ground-state structure at 10 GPa through a fixed-composition structural search. Due to its large unit cell and low symmetry, it can be assumed that the P1 phase, relaxed at the different pressures, approximates quite well the true ground-state of the system; indeed, at all pressures the P1 BaSiH8 phase lies only a few meV above the hull. In practice, analyzing the VCNEB images in Fig. 6, we observe that the transition from the \(Fm\bar{3}m\) to the P1 phase corresponds to the decomposition of BaSiH8 into BaSiH6 + H2, with the expulsion of hydrogen in molecular form.

Fig. 6: VCNEB calculations for BaSiH8 at different pressures.
figure 6

Minimum energy path (MEP) for BaSiH8 at various pressures, calculated using the VCNEB method. Each point represents an individual crystal structure (image) along the path. Particularly relevant intermediate states are shown.

Full size image

The potential barrier, shown in Supplementary Fig. 3, decreases with pressure from 153 meV/atom at 100 GPa to 57 meV/atom at 35 GPa; a sharp transition is visible at 25 GPa, where the kinetic barrier abruptly drops to 9 meV/atom. Although a small barrier survives down to 10 GPa, this sharp decrease is the signature of an impending kinetic instability, i.e., the metastable state will be short-lived and most likely will not be observed in experiments.

Combining the convex hull results in Fig. 1 with the VCNEB analysis, we can argue with confidence that an \(Fm\bar{3}m\) BaSiH8 phase could be synthesized above 100 GPa, and retained down to ~30 GPa, where a clearly visible enthalpy barrier exists. At lower pressures, metastable \(Fm\bar{3}m\) BaSiH8 will decompose, even though (anharmonic) lattice dynamics calculations predict it to be stable. Hence, kinetic stability poses a stricter bound for synthesizability than dynamical stability.

In addition to the promising theoretical synthesizability, both materials appear to be convenient in the experimental setting. As the ambient-pressure convex hulls in Figs. S7 and S9 reveal, both systems feature known stable orthorhombic binary monosilicides, namely BaSi and SrSi43,44, which serve as adequate starting materials having the target composition of Ba/Sr and Si. The monosilicide can be loaded in a diamond anvil cell together with the common solid hydrogen storage medium ammonia borane (H3NBH3) which releases the hydrogen upon heating and forms the refractory compound boron nitride (BN).

Discussion

In summary, our study shows that indeed the superconducting properties of the XYH8 template can be optimized by a suitable choice of the X and Y elements. The two silicides identified in this work represent an improvement compared to LaBH8: for SrSiH8 we predict a dynamical stability pressure pc of 27 GPa, with a Tc of 127 K. Using the figure of merit S introduced by the authors of Ref. 17, this means passing from 1.3 in H3S and LaH10 to 2.2 in LaBH8 to 2.7 in SrSiH8. Even more remarkably, BaSiH8 is predicted to be dynamically stable down to pressures of 3 GPa, with a critical temperature of 71 K, which is substantially higher than all established45 and claimed46 experimental Tc records for conventional superconductors at ambient pressure.

VCNEB calculations demonstrate that the shape of the potential energy landscape of BaSiH8 is favorable for its synthesis in a wide range of pressures. In contrast to the binary La-B system, the Ba-Si system features also the ideal starting material for the synthesis, namely the stable monosilicide BaSi, and a possible path involves synthesis at high pressure (p > 130 GPa) and/or laser heating, and rapid quenching of the resulting phase to lower pressure, i.e., down to ~30 GPa, where the metastable \(Fm\bar{3}m\) crystal structure is protected by a sizable potential barrier. This defines a kinetic threshold pressure, pkin, which is substantially higher than pc.

By defining a concrete method to estimate the synthesizability of a proposed structure, our study sets a new standard for the ab-initio design of new superconductors at high pressures, based on the more rigorous concept of kinetic stability, rather than dynamical stability. The existence of a distinct kinetic stability criterion may also be invoked to explain why many long-standing predictions of high-Tc superconductors have not been realized experimentally47,48,49.

We strongly believe that the proposed method represents a major step forward towards achieving high-Tc conventional superconductivity at room pressure.

Methods

Crystal structure prediction

Crystal structure prediction runs were carried out using evolutionary algorithms as implemented in the USPEX package30,31,50; the underlying total energy and structural relaxation calculations were performed using plane-waves and pseudopotentials as implemented in the Vienna ab-initio Simulation Package VASP51. Further computational details are provided in Supplementary Method 2 of the SM.

Electronic and vibrational properties

All density-functional theory (DFT) and density-functional perturbation theory (DFPT) calculations of electronic and vibrational properties were carried out using the plane-wave pseudoptential code QUANTUM ESPRESSO52, scalar-relativistic optimized norm-conserving Vanderbilt pseudopotentials (ONCV)53, and the PBE-GGA exchange and correlation functional54. Computational details are provided in Supplementary Method 1 of the SM.

Phase transition paths

The phase transition path between the P1 and \(Fm\bar{3}m\)-BaSiH8 was evaluated using the VCNEB method as implemented in USPEX24,30, using variable elastic constants and a variable number of images between the endpoints. The energy and forces were calculated using VASP. Further computational details are provided in Supplementary Method 3 of the SM.

Migdal-Eliashberg theory

The interpolation of the ep matrix elements onto dense wave-vector grids and the subsequent self-consistent solution of the fully anisotropic Migdal-Eliashberg equations were done with EPW41. Based on our previous work on LaBH8, where we calculated the Morel-Anderson pseudopotential μ* from first principles using GW and found consistent values for μ* of about 0.119,25,38, we chose the same value in the current work. Further computational details are provided in Supplementary Method 5 of the SM.

Phonon anharmonicity

Anharmonic corrections to phonon frequencies were obtained by explicitly solving the Schrödinger equation for the APES of every soft mode on a 2 × 2 × 2 wave-vector grid. By calculating and solving 2D APES, we also checked that phonon-phonon interactions for various modes and wave vectors can be neglected in a good first approximation. The interpolation of the real-space force constants obtained on the 2 × 2 × 2 q-grid was performed using the corresponding harmonic support DFPT solution as implemented in the CELLCONSTRUCTOR package55. Further computational details are provided in Supplementary Method 4 of the SM.