In-silico synthesis of lowest-pressure high-$T_\text{c}$ ternary superhydrides

We report the theoretical prediction of two high-performing hydride superconductors BaSiH$_8$ and SrSiH$_8$. They are thermodynamically stable above pressures of $130$ and $174\,\text{GPa}$, respectively, and metastable below that. Employing anharmonic phonon calculations, we determine the minimum pressures of dynamical stability to be around $3\,\text{GPa}$ for BaSiH$_8$ and $27\,\text{GPa}$ for SrSiH$_8$, and using the fully anisotropic Migdal-Eliashberg theory, we predict $T_\text{c}$'s around $71$ and $126\,\text{K}$, respectively. We also introduce a method to estimate the lowest pressure of synthesis, based on the calculation of the enthalpy barriers protecting the BaSiH$_8$ $Fm\bar{3}m$ structure from decomposition at various pressures. This $kinetic$ pressure threshold is sensibly higher than the one based on $dynamic$ stability, but gives a much more rigorous limit for synthesizability.

While it is very exciting to strive for new superconductors with even higher T c 's attaining hot superconductivity [14][15][16], lowering the required stabilization pressures while retaining T c above the boiling point of nitrogen (77 K) is of even greater importance [16][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 larger [6,17]. We have recently predicted, for example, that a ternary hydride with LaBH 8 composition could form an F m3m ternary sodalite-like clathrate (SLC) structure and remain stable down to a critical pressure p c ∼ 35 GPa, with a T c = 126 K [16,19]. Our findings were later confirmed by independent studies on the La-B-H system [21,23]. The predicted p c of LaBH 8 is a factor of four lower than in binary hydrides, where p c '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 F m3m XY H 8 structural template, the discovery of LaBH 8 paved the † lilia.boeri@uniroma1.it * christoph.heil@tugraz.at road to the study of a whole new family of potential high-T c ternary hydrides, where p c and T c may be improved even further. In this work, using first-principles methods for crystal structure prediction and superconductivity, we identify two high-T c alkaline-earth/silicon superhydrides, BaSiH 8 and SrSiH 8 . We predict that both compounds will spontaneously form in the F m3m XY H 8 structure at high pressures (130 and 174 GPa, respectively), and remain dynamically stable down to much lower pressures (p c = 3 and 27 GPa), with superconducting T c 's of 71 and 126 K, respectively.
We also succeed in deriving a rigorous limit for the stability of BaSiH 8 , calculating explicitly the energy barrier protecting the metastable F m3m structure from decomposition (kinetic stability) as a function of pressure, using the Variable-Cell Nudged Elastic Band (VCNEB) method [24]. We find that, indeed, F m3m BaSiH 8 should remain metastable at pressures well below 130 GPa. However, the kinetic critical pressure p kin 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 experiments [11,[25][26][27].

RESULTS
Initial screening: BaSiH 8 and SrSiH 8 were first identified through a high-throughput screening of possible substitutions of La and B by neighbouring elements into the F m3m XY H 8 structure, as shown in Fig. 1a. In this structure - Fig. 1b - cages. This structure permits a very efficient realization of a mechanical precompression mechanism observed in SLC binary hydrides [25,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 p c [19].
We assumed that, considering different combinations of large (X) and small (Y ) atoms, the superconducting properties of LaBH 8 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 F m3m XY H 8 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. Figure 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 = R X + R Y [29], with slopes determined by the total number of valence electrons N e = N X + N Y . Out of the 54 compounds investigated, only three compounds with N e = 6 -green ticks in Fig. 1a -are dynamically stable at 50 GPa: Apart from the already reported LaBH 8 [19], we also found two silicon hydrides, BaSiH 8 and SrSiH 8 . Both compounds exhibit a larger R than LaBH 8 , suggesting that mechanical precompression in both compounds will be more efficient, and hence their p c may be lower. In particular, BaSiH 8 , where R is considerably larger than in the other two compounds, looks very promising in this respect. Indeed, our harmonic phonon calculations yield p c = 5 and 30 GPa for BaSiH 8 and SrSiH 8 , respectively, lower than the p c = 40 GPa in LaBH 8 .

As shown in Supplementary
Convex hulls and phase diagrams: Having established that BaSiH 8 and SrSiH 8 are dynamically stable in the F m3m structure close to ambient pressure, we determined the pressures at which the two compounds may spontaneously crystallize in the F m3m 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 Uspex [30,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 hydrides [32,33]. Further computational details are given in Supplementary Method 2 and Supplementary Note 2 in the SM. An example is illustrated in Fig. 1c In the Ba-Si-H system, the BaSiH 8 composition is thermodynamically stable, i.e., it lies on the hull at 100 GPa in a distorted P 1 phase. The F m3m 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, SrSiH 8 should decompose into SrSiH 6 and SrSiH 12 , which are the stable ternary compositions along with Sr 2 SiH 10 . The corresponding structures are shown in Supplementary Figure 12 of the SM. However, the F m3m 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 BaSiH 8 and SrSiH 8 experience remarkable phonon softening at lower pressures. More specifically, we find soft phonon modes, mainly at Γ (T 2g , E g , and A 1g ) and X (E g , and additionally A 2g in SrSiH 8 ), as shown in Fig. 2. Eventually, imaginary frequencies appear for the E g mode at X below 5 GPa for BaSiH 8 and 40 GPa for SrSiH 8 . (A full set of phonon dispersions at all pressures considered is provided in the SM.) Soft-mode behaviour, associated with strong anharmonicity, has been reported for many hydrogen-rich materials [25,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 BaSiH 8 and SrSiH 8 , 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 BaSiH 8 and SrSiH 8 close to their critical pressures of stability are shown as dashed and full lines in Fig. 2. Anharmonicity causes a considerable hardening of the T 2g (at Γ) and E g (at Γ and X) modes, leading to a decrease of the critical pressure of dynamical stability p c : Taking this hardening into account, SrSiH 8 is stable down to a pressure of 27 GPa and, even more excitingly, BaSiH 8 down to 3 GPa.
Electronic structure and superconductivity: Having determined that both BaSiH 8 and SrSiH 8 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 BaSiH 8 (left) and SrSiH 8 (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 BaSiH 8 (−3.5 meV for SrSiH 8 ), the bands are of purely hydrogen character. This is an important prerequisite for high-T c conventional superconductivity in hydrides, since it can imply a strong electron-phonon (ep) coupling between electronic states to the light hydrogen sublattice [38].
Similar conclusions can be inferred from the electron localization function (ELF) of the two compounds (see Supplementary Figure 15 in the SM), where, in line with what we observed in LaBH 8 [19], 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 hydrides [6,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).
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 BaSiH 8 at 5 GPa pressure and for SrSiH 8 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 code [40,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 BaSiH 8 and SrSiH 8 , respectively.
Further details on the origin of the remarkable T c 's of BaSiH 8 and SrSiH 8 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 α 2 F (ω), 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 BaSiH 8 , for example, we estimate that the phonons related to Γ − T 2g contribute roughly 40 % to the total λ and around X − E g 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 T 2g mode at Γ is a centrosymmetric stretching and squeezing along one of the space diagonals of the H cube, while the E g mode at X is a non-deforming rotation around a face normal of the H cube, with alternating phase in two neighbouring cells.

Metastability:
In all theoretical studies of highpressure 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 values [11,[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 phases [42]. Attempts to quantitatively estimate the barrier height are extremely rare. For F m3m BaSiH 8 , where the dynamical instability pressure p c 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 code [24]. 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 BaSiH 8 at six different pressures: 10, 25, 35, 50, 100, and 200 GPa. As end-members for the VCNEB path, we chose the F m3m phase and a P 1 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 P 1 phase, relaxed at the different pressures, approximates quite well the true ground-state of the system; indeed, at all pressures the P 1 BaSiH 8 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 F m3m to the P 1 phase corresponds to the decomposition of BaSiH 8 into BaSiH 6 + H 2 , with the expulsion of hydrogen in molecular form.
The potential barrier, shown in Supplementary Figure 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 F m3m BaSiH 8 phase could be synthesized above 100 GPa, and retained down to ∼30 GPa, where a clearly visible enthalpy barrier exists. At lower pressures, metastable F m3m BaSiH 8 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 SrSi [43,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 (H 3 NBH 3 ) 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 XY H 8 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 LaBH 8 : for SrSiH 8 we predict a dynamical stability pressure p c of 27 GPa, with a T c of 127 K. Using the figure of merit S introduced by the authors of Ref. [17], this means passing from 1.3 in H 3 S and LaH 10 to 2.2 in LaBH 8 to 2.7 in SrSiH 8 . Even more remarkably, BaSiH 8 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 established [45] and claimed [46] experimental T c records for conventional superconductors at ambient pressure.
VCNEB calculations demonstrate that the shape of the potential energy landscape of BaSiH 8 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 F m3m crystal structure is protected by a sizable potential barrier. This defines a kinetic threshold pressure, p kin , which is substantially higher than p c .
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 predic-tions of high-T c superconductors have not been realized experimentally [47][48][49].
We strongly believe that the proposed method represents a major step forward towards achieving high-T c conventional superconductivity at room pressure.

Crystal structure prediction
Crystal structure prediction runs were carried out using evolutionary algorithms as implemented in the Uspex package [30,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 Vasp [51]. Further computational details are provided in Supplementary Method 2 of the SM.

Electronic and vibrational properties
All density-functional theory (DFT) and densityfunctional perturbation theory (DFPT) calculations of electronic and vibrational properties were carried out using the plane-wave pseudoptential code Quantum Espresso [52], scalar-relativistic optimized normconserving Vanderbilt pseudopotentials (ONCV) [53], and the PBE-GGA exchange and correlation functional [54]. Computational details are provided in Supplementary Method 1 of the SM.

Phase transition paths
The phase transition path between the P 1 and F m3m-BaSiH 8 was evaluated using the VCNEB method as implemented in Uspex [24,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 Epw [41]. Based on our previous work on LaBH 8 , where we calculated the Morel-Anderson pseudopotential µ * from first principles using GW and found consistent values for µ * of about 0.1 [19,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 wavevector 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 package [55]. Further computational details are provided in Supplementary Method 4 of the SM.

DATA AVAILABILITY
The authors confirm that the data supporting the findings of this study are available within the article and its supplementary materials. Further information is available upon request.