Assessing the feasibility of near-ambient conditions superconductivity in the Lu-N-H system

The recent report of near-ambient superconductivity in nitrogen-doped lutetium hydrides (Lu-N-H) has generated a great interest. However, conflicting results have raised doubts regarding superconductivity. Here, we combine high-throughput crystal structure predictions with a fast predictor of the superconducting critical temperature ($T_c$) to shed light on the properties of Lu-N-H at 1 GPa. None of the predicted structures shows the potential to support high-temperature superconductivity and the inclusion of nitrogen favors the appearance of insulating phases. Despite the lack of near-ambient superconductivity, we consider alternative metastable templates and study their $T_c$ and dynamical stability including quantum anharmonic effects. The cubic Lu$_4$H$_{11}$N exhibits a high $T_c$ of 100 K at 20 GPa, a large increase compared to 30 K obtained in its parent LuH$_3$. Interestingly, it has a similar X-ray pattern to the experimentally observed one. The LaH$_{10}$-like LuH$_{10}$ and CaH$_6$-like LuH$_6$ become high-temperature superconductors at 175 GPa and 100 GPa, with $T_c$ of 286 K and 246 K, respectively. Our findings suggest that high-temperature superconductivity is not possible in stable phases at near-ambient pressure, but metastable high-$T_c$ templates exist at moderate and high pressures.


I. INTRODUCTION
that these 1-GPa-phases are improbable to manifest high-temperature superconductivity, as deduced by predicting their T c s with the networking value model [42]. Our findings suggest that, when seeking to develop metallic lutetium hydrides at 1 GPa, doping with nitrogen should be avoided as its large electronegativity removes electrons from hydrogen sites and promotes insulating phases. Moreover, we identify tens of stable metallic phases that exhibit XRD features strongly resembling the experimental XRD, suggesting that many of the predicted structures have a Lu arrangement not far from the fcc lattice and that H or N atoms occupy interstitial sites. Considering that XRD is not capable of distinguishing them, one should approach XRD structural assignments in the Lu-H-N system with care.
As a result of the absence of high-temperature superconductivity at 1 GPa, we study the dynamical stability and superconducting properties of high-symmetry lutetium hydrides with and without nitrogen at higher pressures in crystal structures that favor high T c s. We find that quantum anharmonic effects foster dynamical stability at lower pressures in all cases and strongly impact the phonon spectra. Specifically, cubic Lu 4 H 11 N exhibits a high T c of 100 K at a moderate pressure of 20 GPa. Upon increasing pressure, CaH 6 -like Im3m LuH 6 and LaH 10 -like F m3m LuH 10 are found to maintain high T c s of 246 K and 289 K, respectively, at 100 GPa and 175 GPa.

A. Phase diagram
The phase diagram for the Lu-N-H system at 1 GPa is constructed by the convex hull in The P 2 1 /m LuH 2 N, with structure ID in our database of 2fu_LuH2N_389, is identified as the lowest-enthalpy state among the ternary Lu-H-N compounds at 1 GPa. Dynamical stability in the harmonic approximation of the P 2 1 /m LuH 2 N is examined by performing finite displacement DFT calculations of 2×2×2 supercells. The phonon spectra and phonon density of states (DOS), as depicted in Fig. S1, demonstrate the absence of any imaginary modes, thereby substantiating the dynamic stability of the system. Seeing that the P 2 1 /m LuH 2 N is found to be thermodynamically metastable near the convex hull, it is likely to be synthesized under proper experimental conditions. However, the electronic DOS of the P 2 1 /m LuH 2 N as shown in Fig. S2 exhibits a band gap of ∼ 2 eV, undoubtedly excluding it to be a superconducting phase. Because P 2 1 /m LuH 2 N is the lowest-enthalpy state and is dynamically stable, the phase diagram with respect to stable P3c1 LuH 3 and metastable LuH 2 N is explicitly shown in Fig. 1(b), in which the structures up to H hull of 0.3 eV·atom −1 are displayed. Fig. 1(c) shows the phase diagram with respect to elements, in which the binary hydrides with H hull ≤ 0.8 eV·atom −1 are included. The F m3m phase of LuH 3 , which has been proposed by Dasenbrock-Gammon et al. [26] as the parent compound of the roomtemperature superconductor, is however found to be located above the convex hull by around 82 meV·atom −1 . This enthalpy difference at this pressure is not expected to be overcome by ionic zero-point energy even if quantum anharmonic effects are considered [18].
To identify potential hosts for superconductivity, we have focused our attention on the stability of the 214 metallic phases with H hull ≤ 0.24 eV·atom −1 that have been found in our high-throughput crystal structure predictions. In order to assess the phonon stability of the 214 metallic phases in the harmonic approximation, we have performed high-throughput DFT calculations over 165,000 supercells generated by the finite displacement method (see  To gain more insights into the possible onset of superconductivity among the metallic phases, an understanding beyond H DOS is necessary. However, performing electron-phonon coupling calculations for all these low-symmetry systems is not feasible. As an alternative, Belli et al. have proposed that a physical quantity termed as the networking value (φ), which is the electron localization function value that creates an isosurface spanning throught the whole crystal, can be used to predict easily T c . In fact, φ exhibits a stronger correlation with the actual T c of hydrogen-based superconductors than any other descriptor used so far, and can be used to predict T c with an accuracy of about 60 K with the formula where H f is the hydrogen fraction in the compound [42].
The results of φ and the estimated T c of the dynamically stable phases and all the 214 metallic states irrespective of the stability are included in Table I and Table S1, respectively (an Excel file is also available in Data Availability). Our results indicate that the upper limit of the predicted T c among the 214 metallic states is only 13.94 ± 60 K. These results collectively indicate that the metallic phases predicted from the high-throughput crystal structure prediction are unlikely to host high-temperature superconductivity at 1 GPa.   This particularly short N-H distance corresponds to the strong covalent bonding between N and H, which is evidenced by the electron localization functions of metallic 1fu_LuH2N_366 and insulating 1fu_Lu4H11N_251 belonging to cluster 2 in the lower panel of Fig  This also implies that it is unlikely to find high-temperature superconductivity in the structures from the high-throughput structure screening at 1 GPa.

C. XRD comparison at 1 GPa
Among the binary lutetium hydrides that have been studied experimentally, F m3m LuH 2 has been suggested by Xie et al. [44] and Ming et al. [29] to have the most similar XRD pattern as the one measured by Dasenbrock-Gammon et al. [26]. We thus compare the XRD of LuH 2 with all 638 predicted phases from our high-throughput screening with H hull ≤ 0.24 eV·atom −1 at 1 GPa. In order to compare the XRDs of different structures quantitatively, the similarity between two XRDs is computed according to the correlation function implemented in PyXtal [45]. For simpler comparison, the XRD of F m3m LuH 2 is set as the reference and its simulated XRD is shown in Fig. 3 as a comparison with the experimental XRD reported by Dasenbrock-Gammon et al. [26]. The XRD similarity percentages (%) for the 57 dynamically stable metallic structures are shown in Table. I, while a comprehensive list of the 214 metallic states, irrespective of their dynamical stability, is available in Table S1.  Fig. 3, the dynamically stable binary hydrides with XRD similarity percentages larger than 99% are presented, together with the dynamically stable Cm Lu 4 H 7 N (ID: 1fu_Lu4H7N_11) displaying XRD similarity higher than 90%. In addition, the crystal structures are also shown in the inset of Fig. 3.
Observing the space groups and the crystal structures in Fig. 3, although all six metallic phases show very high XRD similarity, the crystal structure varies from cubic lattice with high symmetry of P m3m to the monoclinic lattice with low symmetry of Cm. The high XRD similarity can even occur in the triclinic structures. For instance, the 2fu_Lu4H11N_570 structure, included in Table S1, shows a large XRD similarity of 93.23% despite its space group is P 1.
The pronounced similarity in XRD can be attributed to the Lu sublattice, which does not differ significantly from the fcc one even if the symmetry reduction is considerable. In compared with the reference state F m3m LuH 2 and the experimental XRD in Ref. [26]. Structure ID, space group, and XRD similarity with respect to F m3m LuH 2 are displayed, along with crystal structure maps.
Figs. S3, S4 and S5 which show the phonon spectra and XRD similarity of the 57 dynamically stable metallic phases, we find that 20 structures exhibiting high XRD similarity greater than 95% share a common feature that their distribution of phonons in frequencies resemble that of F m3m LuH 2 at 1 GPa displayed in Fig. S6a. The specific characteristic is that most phonon branches are separated into two distinct frequency regions, in which some phonons are distributed below 6 THz, which have mainly a Lu character, and other branches are located around 32 THz. By comparing these spectra with the one of LuH 2 (see Fig. S6b), we find that all frequencies around 32 THz in these compounds are associated to vibrations of H atoms around tetrahedral interstitial sites. This means that, despite the notable variations in space groups among the 20 metallic states in Table I In the harmonic approximation, this Lu 4 H 11 N, whose space group is P m3m, is unstable at 1 GPa. The other two high-symmetry binary structures considered are LuH 6 and LuH 10 , which are artificially constructed based on the already known high-temperature hydrogenbased superconductors Im3m CaH 6 [13] and F m3m LaH 10 [18]. It is noted that the crystal structure of Im3m LuH 6 has been theoretically reported in Ref. [46]. At 1 GPa, the phase diagram in Fig. 1c shows that Im3m LuH 6 and F m3m LuH 10 are located around 0.64 and 0.77 eV·atom −1 above the convex hull, implying a highly unstable nature at 1 GPa and 0 K. They are also both dynamically unstable at this pressure at the harmonic level. The zero-point energy is not enough to make any of these structures energetically competitive at 1 GPa.
In order to investigate the impact of quantum anharmonic effects on the dynamical stability of these high-symmetry structures, we relax them within the stochastic self-consistent harmonic approximation (SSCHA) [47][48][49][50] at 300 K at different pressures. This completes the prior study performed on the dynamical stability of the potential parent F m3m LuH 2 and F m3m LuH 3 phases [27]. In order to assess the dynamical stability of these phases we calculate the phonons derived from the Hessian of the SSCHA free energy and check for the presence of imaginary phonon modes [49].
Our SSCHA analysis shows that P m3m Lu 4 H 11 N becomes dynamically stable at around 20 GPa and 300 K (see Fig. 4), which is comparable to the stability range of LuH 3 (6 GPa and 300 K) [27]. The binary compound Im3m LuH 6 has been predicted as a high-temperature superconductor at 100 GPa without the inclusion of ionic quantum anharmonic effects [46]. We fully relax this compound at different pressures within the SSCHA to determine the phase diagram. Although we find it unstable at 100 GPa in the harmonic approximation, it is stabilized with anharmonic and quantum effects at 300 K. The instability in the harmonic approximation is localized at a singular q-point, which in Ref. [46] is also showing significant softening. The q point that shows significant softening in Ref. [46] appears in an 8 × 8 × 8 grid, which could explain why the instability is not present in the previous study that uses instead a 6 × 6 × 6 grid. The calculated electron-phonon coupling constant is huge, which explains the softening of this phonon mode. We find nearly room-temperature superconductivity at 100 GPa with T c of 246 K in this structure, similar to the value reported in Ref. [46], where a T c of 273 K was predicted at 100 GPa at the harmonic level.
The Fm3m LuH 10 structure is the one adopted by the high-temperature superconductor LaH 10 . We find that the free energy Hessian does not show imaginary frequencies above 175 GPa and 300 K indicating a metastable state. We calculate the electron-phonon coupling for this structure and found that the onset of superconductivity happens practically at room temperature, T c = 289 K (∼ 16 • C). In comparison to another polymorph of LuH 10 with a space group P 63/mmc, which can also be stabilized at 200 GPa in the harmonic approximation and has a theoretically reported T c of 134-152 K [51], F m3m LuH 10 phase in our study not only refreshes the record of highest T c of LuH 10 but it can also be stabilized dynamically at a reduced pressure.
In order to confirm the capacity of the networking value model to predict T c s used in the high-throughput calculations, we also estimate T c for these high-symmetry structures with it for the SSCHA structures. We found that the T c of F m3m LuH 3 and P m3m Lu 4

III. CONCLUSIONS
In conclusion, we have performed a comprehensive study by combining a high-throughput structure screening, a rapid estimator of T c , and electron-phonon coupling calculations including quantum anharmonic effects to explore the feasibility of near-ambient superconductivity in the Lu-N-H systems. Our study suggests that the presence of nitrogen favors the formation of insulating states rather than metallic ones, and destabilizes the lutetium hydride systems by significantly shifting them away from the convex hull. As a result, the majority of identified dynamically stable metallic phases in our investigation are binary hydrides, which are not far from the parent F m3m LuH 2 , rather than ternary hydrides. Furthermore, we did not observe high-temperature superconductivity in all of the studied structures at 1 GPa within a reasonable threshold for metastability. We, therefore, propose that in order to have metallic and superconducting states in lutetium hydrides it is better to avoid nitrogen doping. Despite the absence of near-ambient superconductivity, the combined effects of pressure and quantum anharmonicity kindles the hope of high-temperature superconductivity in the Lu-N-H systems by realizing a T c of 100 K at only 20 GPa in cubic Lu 4 NH 11 . This structure is similar to F m3m LuH 3 , but with one out of four H atoms in octahedral sites substituted by a nitrogen atom. Even if this structure is far from being thermodynamically stable at 1 GPa and is dynamically unstable below 20 GPa in our calculations, it is practically indistinguishable from the parent F m3m LuH 2 compound in diffraction experiments. Thus, this structure or variations of it provide, if any, the only possible high-T c structure at low pressures in the Lu-H-N system. At higher pressures, above 100 GPa, CaH 6 -like LuH 6 and LaH 10 -like LuH 10 are superconductors with critical temperatures around room temperature, but with an XRD pattern incompatible with experiments.

A. High-throughput crystal structure prediction
State-of-the-art crystal structure prediction methods, i.e., the evolutionary algorithm implemented in CrySPY [52] and the particle swarm algorithm implemented in CALYPSO [53,54], were combined to predict crystal structures. Crystal structure predictions were performed within first-principles density functional theory (DFT) calculations using the Vienna Ab initio Simulation Package (VASP) [55,56]. The generalized gradient approximation within the parametrization of Perdew, Burke, and Ernzerhof [57] was used with a Hubbard U correction in the Dudarev's form [58] to improve the accuracy of the energies of the Lu f states. An acceptable value of U = 5.5 eV, commonly used to account for the localized f -states of the lanthanide systems [26,59], was used. The plane wave energy cutoff was set to 450 eV during crystal structure predictions. The k-point grid is generated based on the specific structure by Pymatgen [60] with a relatively high grid density of 60 points per Å −3 of reciprocal cell volume. We have screened over 15,000 crystal structures for the Lu-N-H system. To address the electronic properties of the structures within 0.1 eV per atom above the convex hull, the energy cutoff was improved to 550 eV with an improved k-point grid density of 80 points per Å −3 .
To calculate the hydrogen fraction of the total DOS H DOS at the Fermi level, we used Sumo [61] to extract the hydrogen DOS and the total DOS. The TcESTIME code has been used to estimate the superconducting T c based on the networking value model [42]. The XRD simulation and comparison of the XRD similarity were performed by PyXtal [45].
In the phonon calculations for the crystal structures predicted from high-throughput structure screening, the VASP DFT calculations were combined with the supercell and finite displacement methods implemented in Phonopy [62]. Initially, the unit cells from the crystal structure predictions were further optimized with a cutoff energy of 550 eV and a k-point grid density of 130 per Å −3 of reciprocal cell volume until the energy convergence reaches 10 −8 eV and forces of each atom were less than 10 −3 eV·Å −1 . Subsequently, the optimized cells were expanded to supercells for force calculations in DFT. However, considering there were many structures to be examined, it has become infeasible to consider particularly large supercells. It is noted that atomic interactions in most cases were significantly decreased with increased interatomic distances, and thus a cutoff distance of 7.2 Å was considered to be proper in setting up the supercells. Specifically, if the lattice parameter (a, b, or c) of the unit cell was smaller than 7.2 Å, it was expanded twice, otherwise, it remained unchanged.
With this constraint, the 214 metallic phases with H hull ≤ 0.24 eV·atom −1 from the highthroughput crystal structure predictions resulted in more than 165,000 supercells for DFT calculations.

B. Electron-phonon coupling calculations
We relaxed the P m3m Lu 4 H 11 N, Im3m LuH 6 and F m3m LuH 10 using the stochastic self-consistent harmonic approximation method [47][48][49][50]  To calculate superconducting critical temperature for these compounds we performed electron-phonon calculations for the structures obtained through the SSCHA minimization of total free energy using density functional perturbation theory (DFPT) method as implemented in Quantum Espresso [63,64]. Electron-phonon coupling constants were calculated on a 4 × 4 × 4 q point grid for Lu 4  The Eliashberg spectral function was calculated using phonon frequencies obtained from free energy Hessian. The solution of the isotropic Migal-Eliashberg equation was obtained with the cutoff for Matsubara frequencies of 10 times the highest phonon frequency and the reduced Coulomb interaction of µ * = 0.16. The calculation of the electron-phonon coupling for F m3m LuH 3 was done following the parameters previously used in Ref. [27].

V. CODE AVAILABILITY
The high-throughput crystal structural predictions were carried out using the proprietary code VASP [55,56], with the combination of CALYPSO [53,54] and CrySPY [65].       Appendix C: Networking-value-predicted T c of all the metallic phases within 0.24 eV·atom −1 above the convex hull regardless of the dynamical stability.