Potential high-$T_{\mathrm{c}}$ superconductivity in YCeH_x and LaCeH_x under pressure

Because predicting the ternary phase diagram in its entirety is computationally expensive, and our main target is to determine the stability and superconductivity of the hydrogen-rich phase, this research focuses on stable compounds and crystal structures at 100 GPa, 200 GPa, and 300 GPa along a line chosen by YCeH_x and LaCeH_x ($x$ = 2-8, 10, 12, 14, 16, 18, and 20) for the fixed composition search. This approach is currently seen to be effective for alkaline earths and rare earth metal ternary hydrides. Liang et al. (2019a); Shao et al. (2021a); Cataldo et al. (2020) For the stability analysis of YCeH_x and LaCeH_x, the corresponding Y-H, La-H, and Ce-H binary hydrides have been extensively studied. Drozdov et al. (2019); Chen et al. (2021a); Kong et al. (2021) Theoretical calculations show that there is a hydrogen-rich structure $M$H₉ ($M$ = Y, La, and Ce) that can be stabilized below 100 GPa in the three binary systems. In addition, YH₁₀ is capable of exhibiting room temperature superconductivity. Peng et al. (2017) Therefore, this study mainly used these stable binary hydrides as well as single elements as a reference for thermodynamic stability assessment.

The constructed convex hulls of YCeH_x and LaCeH_x at 200 GPa are presented in Fig. 1. Their stable and metastable phases are highlighted by the symbols of circles and red diamond in Fig. 1, respectively. The more detailed results for the stabilites of YCeH_x and LaCeH_x at difference pressures are presented in Figs. S1-S6 in the Supplementary Material (SM). YCeH₅, YCeH₇, YCeH₈, LaCeH₈, and LaCeH₁₈ are found to be stabilized at 100 GPa. When the pressure increases up to 200 GPa, YCeH₅, YCeH₇, and LaCeH₁₈ vanish from the ternary convex, while LaCeH₂₀ and YCeH₁₈ appear in the stable phases of convex hull. In the pressure range of 100-300 GPa, YCeH₂₀ fails to achieve stability. But further increasing the pressure above 300 GPa, $R\bar{3}m$-YCeH₂₀ becomes stable, and its stable pressure interval is very close to that of YH₁₀. Peng et al. (2017) In addition, in the pressure interval of 300-400 GPa, LaCeH₈, LaCeH₂₀, and YCeH₂₀ undergo a phase transition to the more stable $Pmmn$, $P6_{3}/mmc$, and $P4/mmm$ phases, respectively. The ranges of the stabilization pressures of these new predicted phases have been summarized in Fig. 2. To further explore the resistance of these phases to decomposition at high temperatures, we calculated the ternary convex hull at a finite temperature at 200 GPa. We note that the stable phase (YCeH₈, YCeH₁₈, LaCeH₈, and LaCeH₂₀) at this point is also stable at high temperatures, and its specific energy values have been placed in the SM. Therefore, the above-mentioned can be synthesized in the same way as (La,Y)H_x Semenok et al. (2021) by heating LaCe and YCe alloys with NH₃BH₃ at high pressure.

It is interesting to note that most of the above predicted stable phases, except YCeH₅ and YCeH₇, are composed of clathrate structures enclosed by H. For instance, the clathrate structures of stable YCeH₈, LaCeH₈, YCeH₁₈, LaCeH₁₈, and YCeH₂₀ are shown in Fig. 3. Since the hydrogen-rich compounds are more likely to achieve high-$T\rm_{c}$ superconductivity Flores-Livas et al. (2020), the hydrogen-less ones such as YCeH₅ and YCeH₇ will be discussed elsewhere. The stable hydrides of YCeH_x and LaCeH_x with clathrate structures can be roughly divided into three groups.

(i) The first one consists of H₁₈ cages, as found in $P4/mmm$-YCeH₈, $P4/mmm$-LaCeH₈, and $Pmmn$-LaCeH₈, which are derived from their corresponding parent compounds ($I4/mmm$-YH₄, $I4/mmm$-CeH₄, and $I4/mmm$-LaH₄) Peng et al. (2017) with the same cage structure. If half of Ce atoms at the site $2a$ of the $I4/mmm$-CeH₄ phase are orderly replaced by Y or La atoms, the $P4/mmm$ ternary phase would be formed. Peng et al. (2017) The relative enthalpies of these two structures of LaCeH₄ and YCeH₄ are shown in Figs. S1 and S4. We find that the most stable candidate structures prefer the H₁₈ cages. In addition, LaCeH₈ would be transformed into the $Pmmn$ phase when the pressure is above 300 GPa. (ii) The second group consists of H₂₉ cages, as found in the $P\bar{6}m2$-YCeH₁₈ and $Amm2$-LaCeH₁₈ phases. The $P\bar{6}m2$-YCeH₁₈ phase can be derived by replacing half of Ce atoms at the site 2$d$ of $P6_{3}/mmc$-CeH₉ with Y atoms. As for the parent compounds of $P\bar{6}m2$-YCeH₁₈ and $Amm2$-LaCeH₁₈, the previous study by Peng et al. Peng et al. (2017) has shown that only stable $P6_{3}/mmc$-YH₉ and $P6_{3}/mmc$-CeH₉ phases exist in the pressure range of 100-400 GPa, however LaH₉ does not. For the $Amm2$-LaCeH₁₈ phase, only its internal energy plays a superior role to hinder its possible decomposition pathway of 1/6 LaH₄ + 1/6 LaH₁₀ + CeH₉, and thus it cannot be stabilized at higher pressures. Liu et al. (2017); Peng et al. (2017) (iii) The third group consists of the H₃₂ cages, as found in YCeH₂₀ and LaCeH₂₀. The $R\bar{3}m$-YCeH₂₀ and $R\bar{3}m$-LaCeH₂₀ phases exhibit the same spatial structure with LaYH₂₀, Semenok et al. (2021) and all of them are the supercell structures of their parent phases $Fm\bar{3}m$-$A$H₁₀. More specifically speaking, the equivalent primitive unit cell of $R\bar{3}m$-YCeH₂₀ can be obtained by replacing half of Ce atoms with Y atoms in the $1\times 1\times 2$ extension of the primitive unit cell of $Fm\bar{3}m$-CeH₁₀. LaCeH₂₀ undergoes a phase transition from $R\bar{3}m$ to $P6_{3}/mmc$. For YCeH₂₀, it also undergoes a phase transition at 370 GPa and the new high-pressure phase can be derived by replacing half of Ce atoms in the site $4b$ of the parent phase $Fm\bar{3}m$-CeH₁₀ with Y atoms. Considering the chemical similarity of elements, one may choose the stable phases of binary hydrides as the starting structures to accelerate the search for multi-component hydrides.

To further check the lattice dynamic stability of the above predicted YCeH_x and LaCeH_x phases, we have employed the Phonopy Togo and Tanaka (2015) code to initially calculate their phonon band structures, which are presented in Fig. S7 in the SM. Among them, the $Amm2$-LaCeH₁₈ phases cannot maintain dynamic stability in the pressure range of their thermodynamic stability. To make the $R\bar{3}m$-LaCeH₂₀ phase dynamically stable, the pressure needs to be above 200 GPa. In addition, the high pressure applied in the current mainstream experimental studies of metal hydrides falls into a range of 100-300 GPa. As for the phases stabilized above 300 GPa, both $R\bar{3}m$-YCeH₂₀ and $R\bar{3}m$-LaCeH₂₀ have a certain degree of similarity. Therefore, we pay much attention to the $P4/mmm$-YCeH₈, $P\bar{6}m2$-YCeH₁₈, $R\bar{3}m$-YCeH₂₀, $P4/mmm$-LaCeH₈, and $R\bar{3}m$-LaCeH₂₀ phases for their electronic properties and superconductivity.

The electronic band structure, the electronic density of states (eDOS), of $P4/mmm$-YCeH₈, $P\bar{6}m2$-YCeH₁₈, $R\bar{3}m$-YCeH₂₀, $P4/mmm$-LaCeH₈, and $R\bar{3}m$-LaCeH₂₀ at high pressures are presented in Figs. S8-S10 in the SM, respectively. All these compounds in the pressure range of their thermodynamical stabilities possess the metallic features in their energy band structures. For $P4/mmm$-YCeH₈ and $P4/mmm$-LaCeH₈, there are very steep conduction bands along the $Z\rightarrow\Gamma$ direction crossing the Fermi level ($E_{\mathrm{F}}$), resulting in electron pockets near the $\Gamma$ point. For $R\bar{3}m$-YCeH₂₀ and $R\bar{3}m$-LaCeH₂₀, a band inversion can be found around the $\Gamma$ point and near the $E_{\mathrm{F}}$, which is mainly attributed to the H 1s and Ce 4f orbitals, as seen from the atom- and orbital-weighted band structures in Figs. S8 and S9 in the SM, respectively. The eDOS at $E_{\mathrm{F}}$ ($N_{E_{\mathrm{F}}}$) for $P4/mmm$-YCeH₈, $P\bar{6}m2$-YCeH₁₈, $R\bar{3}m$-YCeH₂₀, $P4/mmm$-LaCeH₈, and $R\bar{3}m$-LaCeH₂₀ are 1.72, 1.39, 1.10, 1.34, and 1.26 states/eV/f.u., respectively. As further analyzed by the atom-projected eDOS, the difference in the aforementioned values of $N_{E_{\mathrm{F}}}$ mainly comes from the different contribution of H and Ce atoms to the $N_{E_{\mathrm{F}}}$. The contribution of Ce atoms at $E_{\mathrm{F}}$ is continuously suppressed as the hydrogen content increases, while the opposite trend is observed for H. However, the decrease of Ce is higher than the increase of H, leading to a constant decrease of the total DOS with increasing H content. According to the arguments of Belli et al., Belli et al. (2021) the superconductivity of hydrides maintains a strong positive correlation with the contribution of H in the total DOS at $E_{\mathrm{F}}$ ($H_{\mathrm{DOS}}$). The $H_{\mathrm{DOS}}$ for these five phases are 9.5%, 34.4%, 38.2%, 11.2%, and 34.7%, respectively. This indicates that the H-rich phase has a higher H-derived DOS at $E_{\mathrm{F}}$ although the corresponding total DOS is low, which may suggest that the possibility of high-$T\rm_{c}$ in H-rich phases to some extent. Belli et al. (2021)

From the orbital-weighted electronic band structures of $P\bar{6}m2$-YCeH₁₈, $R\bar{3}m$-YCeH₂₀ and, $R\bar{3}m$-LaCeH₂₀, as shown in Fig. S9, it can be seen that the bands around $E_{\mathrm{F}}$ are mainly contributed by the hybridization between the Ce 4f and H 1s orbitals. We should point out that they were predicted by the GGA-PBE functional. This raises a question how strong the SOC interaction and the correlated effect would be associated with the Ce 4f orbitals in these compounds. So, we further checked the electronic band structures with the SOC interaction, the correction of the on-site Coulomb interaction of Ce 4f orbitals, and their combination, which are presented in Fig. S10 in the SM. We find that both the SOC interaction and the correction of the on-site Coulomb interaction have much weak effect on the bands around $E_{\mathrm{F}}$, which could be neglected. The superconducting transition temperature is largely affected by the electronic properties at the $E_{\mathrm{F}}$. Therefore, the computational methodology employed in the present study would be acceptable to give a reasonable prediction on the superconducting properties of Ce-containing ternary hydrides. The similar treatment with ignoring the SOC interaction and the correlated effect has also been employed in several of the previous theoretical studies Peng et al. (2017); Wang et al. (2021); Li et al. (2019) on the electron-phonon coupling of binary CeH_x. Unfortunately, the combination of the DFT+U method and the density-functional perturbation theory approach in the latest version (v7.1) of QE code does not support the electron-phonon interaction calculations. Because of the above reasons, our superconductivity calculations did not consider the SOC interaction and the strongly-correlated effect.

The phonon dispersions with the mode-resolved electron-phonon coupling (EPC) constant $\lambda_{q\nu}$, atom-projected phonon density of states (PHDOS), and EPC spectra of $P4/mmm$-YCeH₈, $P\bar{6}m2$-YCeH₁₈, $R\bar{3}m$-YCeH₂₀, $P4/mmm$-LaCeH₈, and $R\bar{3}m$-LaCeH₂₀ at high pressures are presented Fig. 4 and Fig. S12 in the SM. For these five phases, the phonons with frequencies above 10 THz are dominated by the H contribution, accounting for more than 70%. The phonons with frequencies below 10 THz are mainly contributed by La, Ce, and Y atoms. Because of the difference in the atomic masses of Y, La, and Ce, the low-frequency range in the phonon spectrum of YCeH_x is slightly wider than that of LaCeH_x. The EPC constants of YCeH₈ at 100 GPa, YCeH₁₈ at 150 GPa, YCeH₂₀ at 300 GPa, LaCeH₈ at 100 GPa, and LaCeH₂₀ at 250 GPa are 0.41, 2.51, 1.02, 0.35, and 0.99, respectively. The EPC constant of YCeH₁₈ is so significant compared to the other four phases. From the mode-resolved EPC constant $\lambda_{q\nu}$, we find that the largest contribution to the total EPC constant of YCeH₁₈ comes from the optical phonon branch in the frequency range of 10 to 20 THz, which contributes about 0.64 (i.e., 25.7% to the total EPC constant). The rare-earth-atom associated vibrational modes with the frequencies of $<$10 THz for YCeH₁₈ also have moderate contribution to the total EPC constant, as compared with the other four phases. In contrast, the contribution of rare-earth atoms to the EPC for YCeH₂₀ with higher H content is only 0.09. When the pressure is increased from 150 GPa to 300 GPa, the phonons dominated by the H-atom-associated vibration in YCeH₁₈ keep hardening, which can been seen from Fig. 4(a) and Fig. S11(c). The increase of pressure leads to a decrease of the EPC constant of YCeH₁₈, so that the EPC constant of YCeH₁₈ at 300 GPa is very close to that of LaCeH₂₀ at 300 GPa. The EPC constants of $A$CeH₈ ($A$ = Y and La) are so small, compared to $A$CeH₁₈ and $A$CeH₂₀, since both the highest phonon frequencies and the H-derived electronic density of states in the former cases are lower than those in the latter cases. These results suggest that once the hydrogen atoms in metal hydrides simultaneously provide sufficient electron density of states at $E_{\mathrm{F}}$ and sufficient phonon density of states in the high-frequency region, a strong EPC strength would be achieved.

Based on the Allen-Dynes-modified McMillan formula Allen and Dynes (1975) with a widely accepted value of the Coulomb pseudopotential $\mu^{*}$ ($\mu^{*}$ = 0.1 used herein), the $T_{c}$ values of $P4/mmm$-YCeH₈ at 100 GPa, $P\bar{6}m2$-YCeH₁₈ at 150 GPa, $R\bar{3}m$-YCeH₂₀ at 300 GPa, $P4/mmm$-LaCeH₈ at 100 GPa, and $P4/mmm$-LaCeH₈ at 250 GPa are predicted to be 4.8, 173.8, 122.1, 1.8, and 116.1 K, respectively. We note that very recently three independent high-pressure experiments Chen et al. (2022a); Bi et al. (2022); Huang et al. (2022) on the La-Ce-H system have been reported, in which the measurement pressures by Chen et al. Chen et al. (2022a) and Bi et al. Bi et al. (2022) were both less than 130 GPa and the observed $T_{\mathrm{c}}$ was around 180 K. One interesting thing is that even though these two research groups Chen et al. (2022a); Bi et al. (2022) used different ratios of the starting materials for synthesis, the final measurement results were very close, suggesting that both of them synthesized probably the same high-temperature superconducting phase. Our calculation results show that at this pressure (130 GPa) the thermodynamically stable phase is LaCeH₁₈, which is not explored further here because the stabilization of this phase may involve anharmonic approximation effect not the focus of this work. The measurement pressure on La-Ce-H by Huang et al. Huang et al. (2022) was in the range of 140-160 GPa. As mentioned above, LaCeH₂₀ can be stabilized in this pressure range. A high-pressure experiment study Chen et al. (2022c) on the Y-Ce-H system has also been reported very recently, in which the observed $T_{c}$ values are between 97-140 K. The predicted $T_{\mathrm{c}}$ value of YCeH₁₈ by either the McM or AD formula shows a good agreement with this experiment data. To a certain extent, these experiment results show strong support to the validity of our theoretical prediction. In short, the high-pressure phases of LaCeH₁₈, LaCeH₂₀, and YCeH₁₈ discovered computationally in this work have all been confirmed in the very recent high-pressure experiments, among which LaCeH₂₀ and YCeH₁₈ are theoretically predicted for the first time.

As presented above, $A$CeH₈, YCeH₁₈ and $A$CeH₂₀ ($A$ = La and Y) have the space groups of $P4/mmm$, $P\bar{6}m2$ and $R\bar{3}m$, respectively, and the rare earth atoms are encased in the H₁₈, H₂₉ and H₃₂ cages, respectively. Herein we performed the COHP and integrated COHP (ICOHP) analyses to further investigate the chemical bonding in these cage structures. The critical bond lengths of some representative H-H, Ce-H, and $A$-H ($A$ - La, Y) atom pairs and their corresponding ICOHP values that can reflect the bonding strength are shown in Figure 5. According to the convention of COHP analysis, the positive and negative values of $-$COHP indicate the bonding and anti-bonding characteristics, respectively. As seen from Fig. 3, the H₁₈ cage in YCeH₈ can be viewed as consisting of eight quadrangles and four hexagons. Accordingly, the H-H bond in the H₁₈ cages can be roughly classified into three groups: i) H-H bonds shared by two neighboring quadrangles, ii) H-H bonds shared by two neighboring hexagons, and iii) H-H bonds shared by a quadrangle and a hexagon.

The H-H bond lengths in the first and second groups are denoted as D1 and D2, respectively, as indicated in Figs. S14(a) and S14(d) in the SM. These two distinct bond lengths D1 and D2 in YCeH₈ and LaCeH₈ at 100 GPa can be sorted in the following descending order: D2(LaCeH₈) $>$ D1(YCeH₈) $\gtrsim$ D1(LaCeH₈) $>$ D2(YCeH₈). Although the H-H bonds with D1 in YCeH₈ and LaCeH₈ have the very close bond lengths, the projected COHPs for them shown in Figs. 5(a) and 5(c) indicate that this H-H bond in LaCeH₈ exhibits a stronger bonding interaction than in YCeH₈ because of the less anti-bonding states appearing just below $E_{\mathrm{F}}$ and a slightly greater absolute value of the corresponding ICOHP in LaCeH₈. However, compared to the $A$-H and Ce-H bonds from their corresponding ICOHP values, the interaction of this H-H bond would be significantly weaker. For YCeH₁₈ and $A$CeH₂₀ with higher H content, the H-H bonds are significantly shrunken because of a higher pressure required to stabilize YCeH₁₈ and $A$CeH₂₀, the ICOHP values of the H-H bonds in YCeH₁₈ and ACeH₂₀ are almost three times than those in $A$CeH₈, and much less (and even no) anti-bonding states appear below $E_{\mathrm{F}}$ for the H-H bonds in YCeH₂₀ (and LaCeH₂₀). From the H₁₈ cages in $A$CeH₈ to the H₂₉ cages in YCeH₁₈ and H₃₂ cages in $A$CeH₂₀, the cage size increases and the number of quadrangles (hexagons) decreases (increases) due to incorporating more H atoms. These results suggest that the H-H bonds within anti-bonding states in $A$CeH₈ could be much stabilized in YCeH₁₈ and $A$CeH₂₀ because they seem ready to be hybridized with the incorporated additional H atoms. In addition, YCeH₂₀ has some remarkable anti-bonding states appearing around -11 eV for the Ce-H atom pair, as seen from Fig. 5(c), which mainly come from the hybridization between the H-1s and Ce-5p_x orbitals, suggesting that the semi-core states (such as 5p orbitals) of Ce in YCeH₂₀ at the high pressure may contribute to the Ce-H bonding.

The crystal orbital bond index (COBI) can be used to quantify the strength of covalent bonds in covalently bonded solid materials. Müller et al. (2021) The integrated-COBI (abbreviated to ICOBI) values for the H-H bonds in these caged structures of $A$CeH₈, YCeH₁₈ and $A$CeH₂₀ are less than 0.15, indicating very weak bonding interactions between H-H. The overall trend is that the ICOBI of H-H bonds increases with the increase of pressure, but the change is less than 0.02, suggesting that the effect of pressure on the H-H bonding is small. The similar trend is also observed from the analysis of ICOHP. One interesting thing is that the bonding strength of the H-H bonds in $R\bar{3}m$-YCeH₂₀ in the pressure range of 300-400 GPa is slightly stronger than the one in $R\bar{3}m$-LaCeH₂₀ according to both the electron localization function (ELF) and ICOHP analyses, although the corresponding ICOBI values in these two structures are quite close to each other and also exhibit a very weak dependence on the pressure.

The five superconducting compounds reported above, which are both thermodynamically and dynamically stable, are mainly of the $A$CeH₈, $A$CeH₁₈ and $A$CeH₂₀ types. The previously reported ScCaH₈ Shi et al. (2021) and YMgH₈ Song et al. (2022a) also have the same structure as $A$CeH₈. For the $A$CeH₈ type, there are three distinct types of H-H bonds, as mentioned above, in which the D1 decreases as the pressure increases, while the D2 exhibits a weaker dependence on the pressure. In additional, the bonding strength of the H-H bonds can also be described as the ELF.

Recently, Belli et al. Belli et al. (2021) have proposed a strong correlation between the $T_{\mathrm{c}}$ value of superconducting hydride and the weak covalent H-H bonds. They also proposed a possible networking value $\phi$, which is defined as the highest value of the ELF that creates an isosurface spanning through the whole crystal in all three Cartesian directions, to describe the nature of the H-H bond in hydrides. As noted by Belli et al., Belli et al. (2021) the precise acquisition of the $\phi$ value is not uniquely defined. From the radial distribution function of H atoms in our studied YCeH_x and LaCeH_x phases as a function of pressure, as presented in Fig. S14 in the SM, we note that within the radial range of $3$ Å there is a main peak representing the most important distribution of H-H bond lengths. In this way, we obtain the $\phi$ value by taking the average of all ELF values at the valley points of the lines plots of ELF along different H-H bonds that correspond to such a main peak. And according to the model proposed by Belli et al., Belli et al. (2021) the $T_{\mathrm{c}}$ of metal hydride can be predicted by the following equation with an accuracy in the range of 60 K:

where $H_{f}$ is the hydrogen fraction and $H_{\mathrm{DOS}}$ is the hydrogen fraction of the total DOS at $E_{\mathrm{F}}$. As presented in Table 1, the network values $\phi$ of our studied phases are all less than 0.5 and they are very close. Consequently, the $T_{\mathrm{c}}$ values predicted by Eq. (9) are more significantly affected by the contribution of $H_{\mathrm{DOS}}$. We also noted that the $T_{\mathrm{c}}$ values predicted by Eq. (9) are still somewhat deviation from those predicted by the AD formula, i.e., Eq. (3). Through the examination on YCeH₁₈ at 150 and 300 GPa, it can be seen that when the pressure increases, the H-H distance decreases, the overall network $\phi$ value exhibits an increasing trend, and the contribution of H to the density of states at $E_{\mathrm{F}}$ also increases. And thus the Eq. (9) based on ELF predicts that the $T_{\mathrm{c}}$ for YCeH₁₈ increases with the increase of pressure, which is opposite to the trend in the prediction of the AD formula. This indicates that one possible limitation in the use of Eq. (9) to predict $T_{\mathrm{c}}$ may be the pressure-dependence of $T_{\mathrm{c}}$. From the predicted results of Peng et al. Peng et al. (2017) for all binary rare-earth metal hydrides, a trend can be seen that the elements containing more f electrons exhibit a suppressing effect on the $T_{\mathrm{c}}$ value. However, the dataset to build Eq. (9) in the work of Belli et al. Belli et al. (2021) contains too few the f-electron-containing cases (i.e., only CeH₉ and PrH₉ included therein). Therefore, the $T_{\mathrm{c}}$ does exhibit a certain relationship with the $\phi$ and $H_{\mathrm{DOS}}$ values, however, the generalization of Eq. (9) to the pressure-dependence and the f-electron-containing systems needs a further improvement.