自1911年荷兰物理学家Onnes首次发现汞的超导现象以来，寻找高温超导体一直是凝聚态物理领域的重要研究课题之一。根据Bardeen-Cooper-Schrieffer （BCS）理论^[1]，超导转变温度T_c与德拜温度成正比，德拜温度与元素的质量成反比，所以轻质元素是潜在的高温超导材料。氢元素作为元素周期表里最轻的元素，理论上最有可能实现高温超导，甚至是室温超导。氢单质在常压下表现为绝缘体，于是Wigner等^[2]根据高压下晶体结构会发生相变这一性质，预言了在极高压力下氢会由非金属相转变为金属相，并且存在高T_c。Dias等^[3]在495 GPa的压力下观测到单质氢出现金属光泽，但该结果存在较大的争议，目前普遍认为在500 GPa以上才有可能实现氢的金属化。

单质氢的金属化需要极高的压力，这无疑是实现高温超导的一大阻碍。为了有效降低氢金属化所需的压力，2004年Ashcroft^[4]提出了“化学预压缩”原理，即在纯氢中引入其他元素“预压缩”氢，得到以氢为主导的材料，使之在相对较低的压力下实现材料的金属化和高温超导。根据这一原理，人们成功预测了一系列氢基高温超导体^[5]，并且其中一些氢化物已在实验室被成功合成。比较有代表性的是共价型H₃S富氢高温超导体，预测其在200 GPa下的T_c为191～204 K^[6–8]。该预测得到了Drozdov等^[9–10]的实验证实，他们观测到在155 GPa下H₃S的T_c为203 K，继而引发了氢化物高温超导体的研究热潮。理论预测表明，LaH₁₀在210 GPa下的T_c为274～286 K^[11–15]，Drozdov等^[16–17]通过实验也证实了约170 GPa压力下LaH₁₀的T_c为250～260 K，进一步增强了人们在氢化物中寻找到室温超导体的信心。实验还合成了其他笼状结构的氢化物，如CaH₆、YH₆和YH₉，在150～300 GPa的压力范围内，其T_c超过200 K^{[12, 18–20]}。此外，还预测了层状类五角石墨烯结构HfH₁₀在250 GPa压力下的T_c约为230 K^[21]。

目前，二元氢化物的研究几乎覆盖了整个元素周期表。三元氢化物由于其化学成分的多样性，具有更丰富的性质和更广阔的探索空间。在之前的研究中，通过将金属元素引入富含H₂分子的二元氢化物，成功预测了Li₂MgH₁₆在250 GPa下的T_c为473 K，即Li₂MgH₁₆为“热”超导体，是现有理论预测T_c的最高记录^[22–23]。此外，通过在二元氢化物中引入第3种元素还可以进一步降低稳定压强。理论预测显示，LaBeH₈的动力学稳定压力可以降至20 GPa，且T_c越过了液氮温区（约190 K） ^[24]。最近，LaBeH₈已在实验中被成功合成^[25]。理论预测还发现，三元氢化物KB₂H₈的动力学稳定压力可降至12 GPa^[26]。因此，三元氢化物中极有可能存在低压下稳定并且具备高温超导性质的超导体。

相比于二元氢化物，三元氢化物由于多了一种元素参与反应，合成难度增加。目前已合成的三元氢化物，如(La,Ce)H₉、(La,Al)H₁₀以及LaBeH₈，均采取了合金掺入H原子的方法^{[25, 27–28]}。H的掺入使得原本没有超导性质的合金转变成了超导体。因此，将H原子掺入现有的化合物，极有可能使其转变为超导相或者提高其T_c。

近期，Zhang等^[29]在Y-Si体系的研究中预测了电子化合物Y₃Si在30和50 GPa下的T_c分别为11.2和14.5 K。除了这一具有代表性的结构外，另外还发现了YSi、YSi₂、YSi₃、Y₅Si₃、Y₂Si 5种稳定的化合物，但它们几乎不具备超导性质。本研究以Y-Si化合物为基本载体，通过掺入氢原子，利用氢的低原子质量和高振动频率提高T_c或使其转变为超导体，采用第一性原理对Y-Si-H体系进行系统的结构预测、电子性质分析和T_c计算。

1.   计算细节

采用基于密度泛函理论的随机结构搜索软件Airss^[30–31]对三元氢化物体系Y_xSi_yH_z （x=1, 2；y=1, 2； z=1, ···, 20）在150 GPa下进行晶体结构搜索，使用第一性原理计算软件包CASTEP ^[32]对预测结构进行结构优化，在同一配比下生成的结构个数不低于300个。选择广义梯度近似的Perdew-Burke-Ernzerhof ^[33–34]作为交换关联函数。首先，开展结构的低精度搜索。在低精度搜索过程中选用超软赝势^[35]，第一布里渊区的K点网格尺寸选为 2π×0.07 Å⁻¹，截断能为300 eV，以达到快速搜索的目的。然后，采用MAGUS软件包^[36–37]进行结构搜索，以相互验证搜索结构的准确性。结果显示，搜索结果与Airss搜索结果基本一致。接着，对配比焓值最低的10种结构进行高精度优化。优化过程中，截断能选取800 eV，第一布里渊区的K点网格尺寸选为2π×0.03 Å⁻¹。最后，采用VASP进行电子结构计算，赝势采用PAW赝势^[38]，K点网格尺寸选为2π×0.02 Å⁻¹，平面波截断能选为800 eV，能量收敛值取为1×10⁻⁷ eV。此外，为了获得压力下的平衡结构，对内部原子位置进行了优化，原子残余力阈值设为5×10⁻³ eV/ Å。对于电子-声子耦合特性和超导性质的计算，采用基于密度泛函微扰理论的Quantum-ESPRESSO软件包^[39–40]，选用超软赝势^[35]，截断能选为80 Ry （约1100 eV），布里渊区的q点网格数量选为5×5×5，K点网格数量选为20×20×20，以保证计算的准确性。

2.   结果与讨论

2.1   结构预测

首先，采用随机结构搜索方法在150 GPa下对Y-Si-H体系进行变组分结构预测。为进一步研究其热力学稳定性，计算了已报道的单质和稳定的二元化合物在高压下的焓值^[29]，构建热力学凸包图，如图1所示。图1中，红色实心点代表凸包上的热力学稳定结构，蓝色空心点代表偏离凸包图15 meV/atom以内的亚稳结构。通常来说，距凸包图50 meV/atom以内的化合物被认为是半稳定化合物，在特定的温度下有可能被合成并稳定存在^[41–42]。在150 GPa下，YSiH₇、YSiH₉、YSi₂H₁₂和YSiH₁₈均落在凸包图上，为热力学稳定化合物，不会分解为单质或者其他化合物。YSi₂H₁₃、YSi₂H₁₄和Y₂SiH₁₇的焓值分别偏离凸包图4、12和6 meV/atom，是热力学亚稳化合物。

图  1  Y-Si-H体系在150 GPa压力下的热力学凸包图（红色实心点表示落在凸包上的热力学稳定的化合物，蓝色空心点表示偏离凸包15 meV/atom以内的亚稳化合物）

Figure  1.  Thermodynamic convex hull diagram of Y-Si-H system under pressure of 150 GPa (The red solid points are represented as thermodynamically stable structures on the convex hull, and the blue hollow points are represented as metastable structures deviating from the convex hull within 15 meV/atom.)

下载: 全尺寸图片 幻灯片

2.2   晶体结构

图2为Y-Si-H体系在150 GPa压力下的热力学稳定或者亚稳的晶体结构，其中，图2(a)、图2(b)、图2(e)、图2(f)分别对应热力学稳定的YSiH₇、YSiH₉、YSi₂H₁₂和YSiH₁₈，图2(c)、图2(d)、图2(g)分别对应亚稳的YSi₂H₁₃、Y₂SiH₁₇和YSi₂H₁₄。这些结构有一个共同的特点，即Si原子附近均存在大量的H原子，Si的配位数均为6以上，并且大部分结构中存在H₂分子单元。前人总结了高压下高T_c富氢化合物的一般规律^{[12, 21]}：(1)具有高对称性结构，(2)不存在H₂或H₃分子单元，(3)费米面处H的电子态密度占比高，(4)费米面处电子与高频声子强耦合。在YSiH₁₈中，每个Si原子连接6个H原子，6个H原子构成八面体包裹Si，所有Si―H键等价，键长为1.46 Å，但因为YSiH₁₈中除Si―H键外还存在大量的H₂分子单元，所以其超导温度不会太高。相反，对于YSi₂H₁₂，每个H原子都与Si成键，不存在H₂分子单元，因此可能具有较高的T_c。

图  2  预测的热力学稳定或亚稳的三元Y-Si-H化合物的晶体结构

Figure  2.  Predicted crystal structure of thermodynamically stable or metastable ternary Y-Si-H compounds

下载: 全尺寸图片 幻灯片

2.3   合成路径分析

计算了Y_xSi_yH_z的生成焓，包括稳定的化合物YSiH₇、YSiH₉、YSi₂H₁₂和YSiH₁₈以及亚稳化合物YSi₂H₁₃、YSi₂H₁₄和Y₂SiH₁₇，以期为实验合成提供有效的指导^[43–44]。以YH₃、Y金属、Si单质和H单质作为反应的前驱物，构建了两条不同的反应路径，式(1)～式(2)为两条反应路径配平后的生成焓方程，${H}^{{f}_{1}}$和${H}^{{f}_{2}}$分别代表路径1和路径2的生成焓。

图3和图4展示了Y-Si-H体系化合物的生成焓与压力的关系。生成焓小于零表明目标化合物可能通过该反应路径合成，生成焓越低，合成的可能性就越大^[44]。图3和图4显示，在不同的目标化合物中，${H}^{{f}_{1}}$和${H}^{{f}_{2}}$均随着压力的增加而逐渐减小，${H}^{{f}_{2}}$始终小于${H}^{{f}_{1}}$，且在绝大多数压力范围内${H}^{{f}_{2}}$和${H}^{{f}_{1}}$为负值。对于稳定的化合物YSiH₇、YSiH₉、YSi₂H₁₂和YSiH₁₈，其分别在50、95、80和80 GPa以上时${H}^{{f}_{1}}$小于零；对于亚稳化合物YSi₂H₁₃、YSi₂H₁₄和Y₂SiH₁₇，其分别在150、90和100 GPa以上时${H}^{{f}_{1}}$小于零。${H}^{{f}_{1}}$和${H}^{{f}_{2}}$均为负值表明目标化合物可能同时通过这两条反应路径合成，此时合成目标化合物的可能性大幅度增加。

图  3  计算的YSiH₇、YSiH₉、YSi₂H₁₂和YSiH₁₈形成焓随压力的变化关系

Figure  3.  Calculated formation enthalpies of YSiH₇, YSiH₉, YSi₂H₁₂ and YSiH₁₈ as a function of pressure

下载: 全尺寸图片 幻灯片

图  4  计算的YSi₂H₁₃、YSi₂H₁₄和Y₂SiH₁₇的形成焓随压力的变化

Figure  4.  Calculated formation enthalpies of YSi₂H₁₃, YSi₂H₁₄ and Y₂SiH₁₇ as a function of pressure

下载: 全尺寸图片 幻灯片

2.4   电子性质

为了进一步了解Y-Si-H体系结构的电子性质，研究其在150 GPa压力下的能带结构和分波态密度（partial density of states, PDOS），结果如图5和图6所示。图6(a)显示，YSiH₁₈的能带穿越费米能级，但是在费米能级附近仅有少量电子，可能不利于提高T_c。图6(b)显示，YSiH₇的能带在费米能级附近存在一个3 eV的能隙，表明YSiH₇是绝缘体，在此条件下不具备导电性质。图6(c)显示，YSi₂H₁₃虽有能带穿过费米能级，但是从PDOS可以看出，费米能级处基本无电子占据，在Γ点附近存在空穴型费米面口袋，口袋顶部接近导带的电子口袋底部，表明YSi₂H₁₃具有半金属特性。因此，对YSiH₇和YSi₂H₁₃不做后续研究。其余结构的能带图均显示有能带越过费米能级。图5(a)显示，YSiH₉的能带沿Γ-A方向存在一个电子口袋，口袋底部靠近费米能级，其他区域的能带则存在较宽的能隙。对于YSi₂H₁₂（见图5(b)），沿着M-D方向能带在费米能级附近发生交叉，并且从态密度图上看，费米能级附近存在电子占据，表明YSi₂H₁₂具有较好的金属性质。此外，YSi₂H₁₂中H原子在费米能级处占据最大的电子态密度（大于0.7），由于H原子在费米能级处占据高的电子态密度对于提高T_c十分有利^[21]，因此YSi₂H₁₂可能是潜在的高温超导体。YSi₂H₁₄的能带结构沿Γ-Z-R方向（见图5(c)）和Y₂SiH₁₇的能带结构沿Y-M-A方向（见图5(d)）均在费米能级附近出现平带，这种电子结构提高了费米能级附近的电子占据，有利于提高T_c。另外，值得注意的是，YSi₂H₁₄和Y₂SiH₁₇结构中Y原子的d轨道的电子态密度在费米面处占主导（约占0.59和0.61）。最近的研究表明，过渡金属的d轨道电子对超导性质的形成也有帮助^[45–46]。

图  5  YSiH₉、YSi₂H₁₂、YSi₂H₁₄、Y₂SiH₁₇在150 GPa压力下的能带结构和分波态密度

Figure  5.  Band structure and discrete electron state density of (a) YSiH₉, (b) YSi₂H₁₂, (c) YSi₂H₁₄, (d) Y₂SiH₁₇ under pressure of 150 GPa

下载: 全尺寸图片 幻灯片

图  6  YSiH₁₈、YSiH₇、YSi₂H₁₃在150 GPa压力下的能带结构和分波态密度

Figure  6.  Band structure and discrete electron state density of YSiH₁₈, YSiH₇, YSi₂H₁₃ under pressure of 150 GPa

下载: 全尺寸图片 幻灯片

2.5   超导性质

Allen等^[47]对Mc Millan方程进行了简化，当$\lambda$ << 1.5时，超导转变温度为

当$\lambda$ > 1.5时

${\overline{\omega }}_{2}$、${\omega }_{\mathrm{l}\mathrm{o}\mathrm{g}}$和$\lambda$可以分别表示为

式中：${\omega }_{\mathrm{l}\mathrm{o}\mathrm{g}}$ 为声子频率的对数平均值，$\lambda$为电声耦合常数，μ^*为有效的库伦作用势，${\overline{\omega }}_{2}$ 为声子频率平方的平均值。

针对上述具有金属性质的结构，计算其声子谱、谱函数和电声耦合常数$\lambda$。对于氢化物，μ^*一般取 0.10～0.13，本研究中μ^*= 0.13，计算得到的T_c列于表1。表1显示，YSi₂H₁₂在100 GPa下的T_c为43.5 K，超过了麦克米兰极限，是该体系超导转变温度的最高值，该超导转变温度与费米能级处H原子的电子占主导有很大的关系。另外，费米能级处Y原子的d轨道电子占主导地位的YSi₂H₁₄和Y₂SiH₁₇也拥有较高的T_c。Y₂SiH₁₇在100 GPa下的T_c可达35.8 K。YSi₂H₁₄的动力学稳定压力可降至40 GPa，并且T_c高达23.8 K，已接近同等压力下Y-Si二元化合物中T_c最高的Y₃Si的2倍（Y₃Si的T_c约为14 K）。由此可见，H的掺入可以有效地提高Y-Si体系的超导转变温度。对于YSiH₁₈，H原子主要以H₂分子的形式出现，致使费米能级处仅有少量氢电子占据，其T_c值仅为1.9 K。

表  1  Y-Si-H体系各结构在不同压力下的H原子的s轨道和Y原子的d轨道在费米面处的态密度（${\mathit{N}}_{\rm{Ef}}$）、平均声子频率的对数（${\mathit{\omega }}_{\rm{log}}$）、电声耦合参数（$\mathit{\lambda }$）和超导转温度（$T_{\mathrm{c}}$）

Table  1.  Density of state (${\mathit{N}}_{\rm{Ef}}$) of the s orbital H and the d orbital Y at the Fermi surface, logarithm of average phonon frequency ${\mathit{\omega }}_{\rm{log}}$, electroacoustic coupling parameter $\mathit{\lambda }$, and superconducting transition temperature $T_{\mathrm{c}}$for each structure of Y-Si-H system under different pressures

Compound	Phase	Pressure/GPa	NEf /(states∙eV−1)		${\mathit{\omega }}_{{\mathrm{{l}{o}{g}}}}$/K	$\mathit{\lambda }$	Tc/K
			s orbital of H	d orbital of Y
YSiH9	P21/m	150	0.15	0.50	951	0.60	15.3
YSiH18	P312	150	0.11	0.11	605	0.43	1.9
YSi2H12	C2/m	150	0.28	0.06	1117	0.68	28.4
YSi2H12	C2/m	100	0.32	0.07	954	0.87	43.5
YSi2H14	C2/m	150	0.10	0.31	842	0.55	9.7
YSi2H14	C2/m	100	0.09	0.33	776	0.58	11.5
YSi2H14	C2/m	40	0.10	0.44	419	0.97	23.8
Y2SiH17	P$\overline{\mathrm{4}}$m2	150	0.13	0.28	873	0.78	31.6
Y2SiH17	P$\overline{\mathrm{4}}$m2	100	0.13	0.30	694	0.92	35.8

下载: 导出CSV 

| 显示表格

影响T_c的因素很多，其中电声耦合常数的大小直接影响T_c。一般来说，电声耦合常数越大，T_c也就越高。为了进一步确定诱导YSiH₉、YSi₂H₁₂、YSi₂H₁₄和Y₂SiH₁₇超导的原因，计算了它们在各自最低动力学稳定压力下的声子谱、声子态密度、谱函数和电声耦合常数，如图7所示。图7(a)、图7(c)和图7(d)显示，YSiH₉、YSi₂H₁₄和Y₂SiH₁₇在高频区域（大于3000 cm⁻¹）存在孤立平坦的声子谱，对应H₂分子的拉伸振动模式，由此证明了H₂分子单元的存在。YSiH₉的声子谱显示，第11支声子频率（声子频率范围500～1000 cm⁻¹）在Γ点和B点附近明显降低，出现了声子软化现象，并且在该声子频率范围内$\lambda$快速增加。另外，高频H₂分子的拉伸振动对$\lambda$基本没有贡献。YSi₂H₁₄和Y₂SiH₁₇的$\lambda$在低频区分别增加了55.7%和46.2%，从声子态密度（phonon density of states，PHDOS）中可知，在该频率区域中分布着Y原子与H原子的振动模式，除H原子外，Y原子也参与到电声耦合作用当中，因此Y原子对超导也起着重要作用。考虑YSi₂H₁₄和Y₂SiH₁₇的费米能级处Y的d轨道电子占据主导作用，推测在费米能级附近Y的d轨道电子占据可能会诱发高的超导转变温度。对于YSi₂H₁₂（见图5(b)），$\lambda$在声子振动的中高频区域快速增加，对应的频率主要来自H原子的振动，因此，YSi₂H₁₂的$\lambda$主要源于中高频H原子的贡献。

图  7  YSiH₉、YSi₂H₁₂、YSi₂H₁₄、Y₂SiH₁₇分别在150、150、40和100 GPa下的声子谱、声子态密度、谱函数${\alpha }^{2}F(\omega)$和电声耦合常数$\lambda$

Figure  7.  Phonon spectra, phonon state density, spectral function ${\alpha }^{2}F(\omega)$, electroacoustic coupling and $\lambda$ at 150, 150, 40 and 100 GPa of (a) YSiH₉, (b) YSi₂H₁₂, (c) YSi₂H₁₄, (d) Y₂SiH₁₇, respectively

下载: 全尺寸图片 幻灯片

2.6   稳定性分析

化合物的性质与组分之间存在密不可分的关系^[48]，对比YSiH₉、Y₂SiH₁₇与YSi₂H₁₂、YSi₂H₁₄的最低动力学稳定压力发现，当Si原子的含量增加时，其动力学稳定压力随之降低，尤其是YSi₂H₁₄，它在40 GPa下仍保持动力学稳定，因此，可以认为Si原子含量的增加可以有效降低动力学稳定所需的压力。

通过计算预测化合物的弹性常数${C}_{ij}$，检验其在最低动力学稳定压力下的弹性稳定性。经计算，YSiH₉在150 GPa下的弹性常数矩阵为

Y₂SiH₁₇在100 GPa下的弹性常数矩阵为

YSi₂H₁₂在100 GPa下的弹性常数矩阵为

YSi₂H₁₄在40 GPa下的弹性常数矩阵为

波恩稳定性判据为^[49]

劳厄群为4/mmm的Y₂SiH₁₇属于正交类晶体结构，拥有6个独立的${C}_{ij}$。从100 GPa下Y₂SiH₁₇的弹性矩阵可以看出，${C}_{ij}$满足上述判据，表明Y₂SiH₁₇在100 GPa下是力学稳定的。YSiH₉、YSi₂H₁₄和YSi₂H₁₂均为单斜晶系，拥有13个独立的弹性常数，研究低对称性晶体时，通常采用检验弹性矩阵的所有特征值均为正来验证弹性稳定性^[49]。YSiH₉、YSi₂H₁₄和YSi₂H₁₂的弹性矩阵满足特征值均为正的判定标准，因此，所预测的化合物在其最低动力学稳定压力下仍保持着力学稳定。

3.   结　论

采用第一性原理计算方法，在150 GPa下对Y-Si-H体系进行了多种化学配比的结构搜索，确定了热力学稳定的YSiH₇、YSiH₉、YSi₂H₁₂和YSiH₁₈结构，以及亚稳的YSi₂H₁₃、YSi₂H₁₄和Y₂SiH₁₇结构。电子性质计算表明，除YSiH₇为绝缘体、YSi₂H₁₃为半金属外，其他氢化物均具有金属性，是潜在的超导体。电声耦合计算表明，在相同压力下，相对于Y-Si二元化合物中具备最高超导转变温度的Y₃Si，YSi₂H₁₄的超导转变温度约为Y₃Si的2倍，因此，H原子的介入有利于提升Y-Si体系的超导转变温度。对YSi₂H₁₄和Y₂SiH₁₇的超导性质研究表明，Y原子的d轨道电子在费米能级处占据主导和Y原子的低频振动强化电子-声子耦合均对超导温度的提升起着重要作用。