Bulk Superconductivity in the Dirac Semimetal TlSb
Yuxing Zhou, Bin Li, Zhefeng Lou, Huancheng Chen, Qin Chen, Binjie Xu, Chunxiang Wu, Jianhua Du, Jinhu Yang, Hangdong Wang, Minghu Fang
BBulk Superconductivity in the Dirac Semimetal TlSb
Yuxing Zhou, Bin Li, Zhefeng Lou, Huancheng Chen, Qin Chen, Binjie Xu, Chunxiang Wu, Jianhua Du, Jinhu Yang, Hangdong Wang, and Minghu Fang
1, 5, ∗ Department of Physics, Zhejiang University, Hangzhou , China New Energy Technology Engineering Laboratory of Jiangsu Province and School of Science,Nanjing University of Posts and Telecommunications, Nanjing , China Department of Applied Physics, China Jiliang University, Hangzhou , China Department of Physics, Hangzhou Normal University, Hangzhou , China Collaborative Innovation Center of Advanced Microstructures, Nanjing University, Nanjing , China (Dated: December 8, 2020)A feasible strategy to realize the Majorana fermions is searching for a simple compound withboth bulk superconductivity and Dirac surface states. In this paper, we performed calculations ofelectronic band structure, the Fermi surface and surface states, as well as measured the resistivity,magnetization, specific heat for TlSb compound with a CsCl-type structure. The band structurecalculations show that TlSb is a Dirac semimetal when spin-orbit coupling is taken into account.Meanwhile, we first found that TlSb is a type-II superconductor with T c = 4.38 K, H c (0) = 148Oe, H c (0) = 1.12 T and κ GL = 10.6, and confirmed it to be a moderately coupled s -wave supercon-ductor. Although we can not determine which bands near the Fermi level E F to be responsible forsuperconductivity, its coexistence with the topological surface states implies that TlSb compoundmay be a simple material platform to realize the fault-tolerant quantum computations. I. INTRODUCTION
Topological superconductors host Majorana fermionsdescribed by a real wave function, providing protectionfor quantum computations [1]. So, realizing topologicalsuperconductivity (TSC) has became one of the most in-teresting topics in the condensed matter physics in thepast decades. According to the discussion in Ref. [2],there are intrinsic and artificial engineered topologicalsuperconductors. For the intrinsic, the topological non-trivial gap function naturally shows up. Sr RuO [3–6]is the first proposed topological superconductor althoughthe existence of chiral p -wave superconductivity (SC) isstill under debate. Cu x Bi Se [7], as the first material toshow SC ( T c ∼ x Bi Se . The nematic SC discovered inCu x Bi Se [9, 10] was also observed in the similar su-perconductors derived from Bi Se , such as in Sr x Bi Se [11], Nb x Cu Se [12]. Sn − x In x Te [13] is another super-conductor upon doping charge carriers into a topologicalcrystalline insulator. In the cleanest sample ( x ∼ T c (1.2 K), a pronounced zero-bias con-ductance peak (ZBCP) similar to that in Cu x Bi Se hasbeen observed by point contact spectroscopy [13]. An-other is the artificial engineered TSC in hybrid struc-tures. According to the idea proposed by Fu and Kane[14], if s -wave pairing is imposed on the topological sur- face states of a three dimensional (3D) TI through super-conducting proximity effect, the resulting superconduct-ing state should be a 2D p -wave SC harboring a Majoranazero mode in the vortex core. Experimentally, proximity-induced SC on the surface of 3D TIs has been studied bymany groups[15–21]. The observation [22] of 4 π -periodJosephson supercurrent in 3D HgTe TI is encouraging,although, it is difficult to elucidate the topological na-ture of the induced 2D SC.However, there are controversies about TSC emergingand the observed ZBCP being a Majorana zero-energymode (MZM) in the doped topological material. To re-alize TSC through a superconducting proximity effect inhybrid structures has many engineering challenges. Afeasible way is to realize TSC in a simple compound, inwhich both the topological surface state and the bulk SCcoexist, thus Majorana fermions emerge at the edge ofsuperconductor. Recently, the observations of MZM inthe core of vortex [23, 24], at the end of the atomic defectline [25], and near the Bi islands [26] and the interstitialFe atoms [27] in the simplest Fe-based superconductorFe y Te . Se . ( T c = 14 K) [28] motivate us to searchfor the similar system. TlSb crystallizes in a cubic CsClstructure with space group P m ¯3 m (No. 221) (as shownin the inset of Fig. 1), from this structure a large num-ber of topological semimetal/metals (TMs) were designed[29], ranging from triple nodal points, type-I nodal lines,and critical type nodal lines to hybrid nodal lines. Forexample, CaTe is a typical type-I nodal line and DiracTM [30]; YIr is a typical triple-nodal-point TM [29], YMgpossesses multiple types of band crossing [29]. Therefore,we tried to grow TlSb crystals for studying its topolog-ical natures and SC, unfortunately, only polycrystallineTlSb samples were obtained. a r X i v : . [ c ond - m a t . s up r- c on ] D ec Tl Sb
Figure 1. (a) Polycrystalline XRD pattern with its refinementprofile at room temperature of TlSb. The inset shows theschematic structure of TlSb, thallium atoms are in gray whilethe antimony atoms are brown.
In this paper, we performed calculations of the elec-tronic band structure, the Fermi surface and the surfacestates on (001) plane, as well as measured resistivity,magnetization and specific heat for the polycrystallineTlSb sample. The band structure calculations show TlSbis a Dirac semimetal with 4-fold degenerate nodes nearΓ and R points. It is also found that TlSb is a type-II superconductor with the superconducting transitiontemperature T c = 4.38 K, the lower critical field H c (0)= 148 Oe, and the upper critical field H c (0) = 1.12 Tand the Ginzburg-Landau (GL) parameter κ GL = 10.6.The obtained specific heat jump, ∆ C el / γ n T c ∼ s -wave superconducting symmetry.These results indicate that both s -wave SC and surfacestates coexist in TlSb, whether the Majorana fermionsemerge or not on the edges is needed to confirm in thefuture. II. EXPERIMENTAL
Polycrystalline TlSb samples were synthesized by aperitectic reaction method. The mixture of stoichiomet-ric high purity Tl (99.99%) chunk and Sb (99.999%) pow-der was placed in an alumina crucible, sealed in an evac-uated quartz tube and heated at 450 ◦ C for 20 hrs, thendecreases to 195 ◦ C waiting for the mixture melting com-pletely. To avoid the decomposition of obtain TlSb phaseat 191 ◦ C [31], the quartz tube was quenched to roomtemperature at 195 ◦ C, however, the obtained TlSb sam-ples always contain a small amount of unreacted Sb im-purities due to the precipitation of Sb before the peritec-tic reaction. The obtained TlSb alloy is easily to cut forthe subsequent structure characterizations and property measurements. Polycrystalline x-ray diffraction (XRD)was carried out on a PANalytical diffractometer equipedwith CuK α radiation. The TlSb XRD pattern is shownin Fig. 1, in which the main peaks can be fitted by theCsCl-type structure with space group P m ¯3 m . The lat-tice parameters a = b = c = 3.86(5) ˚A were obtained bythe Rietveld refinement by using general structure anal-ysis system (GSAS) [32]. A rectangular bar of the sam-ple was cut for the magnetization and resistivity mesure-ments, which were performed on a magnetic propertymeasurement system (Quantum Design, MPMS - 7 T)and a physical property measurement system (QuantumDesign, PPMS - 9 T), respectively. The band structurewas calculated by using density function theory (DFT)with the WIEN2k package [33]. Generalized gradient ap-proximation (GGA) of Perdew-Burke-Ernzerhof (PBE)[34] was employed for the exchange correlation potentialcalculations. A cutoff energy of 520 eV and a 13 × × k -point mesh were used to perform the bulk calcula-tions. The Fermi surface (FS) was performed with Wan-nierTools [35] package which is based on the maximallylocalized Wannier function tight-binding model [36–38]constructed by using the Wannier90 [39] package. III. RESULTS AND DISCUSSIONS
We first discuss the electronic band structure withoutconsidering spin-orbit coupling (SOC). As shown in Fig.2(a), both conduction (red) and valence (blue) band crossthe Fermi level E F . At the high symmetry Γ and Rpoints, there are three bands crossing, with a threefolddegenerate at 0.8 eV and 2 eV below E F , respectively.However, when SOC is taken into account [see Fig. 2(b)],gaps open and leaves two twofold degenerate bands atboth high symmetry points. Since both time-reversaland inversion symmetries are present, no spin-splittingoccurs and then the twofold degenerated bands come to-gether to fourfold degenerate points, indicating TlSb isa Dirac semimetal. We also calculated the density ofstate (DOS), as shown in the right panel of Fig. 2(b),the DOS at E F is mainly contributed by Sb orbits. Tofurther clarify the band structure of TlSb, we calculatedits 3D bulk FS of the first Brillouin zone (BZ) as shownin Fig. 2(d), exhibiting very complex 3D characteristics.Figure 2(e) presents the FS on the k z = 0 plane, which isthe cross section passing the Γ point of the 3D FS. Figure2(f) displays the energy dispersion in the k x - k y plane, inwhich the Dirac dispersion is clearly seen, demonstratingfurther that TlSb is a Dirac semimetal. Then we calcu-lated the surface states on (001) plane by using a surfaceGreen’s function method [40]. As shown in Fig. 3(a), theprojected Dirac points are hidden in the continuous bulkstates, the surface states are shown as the red curves.The (001) surface energy contour is shown in Fig. 3(b),(c) and (d) with E = E F , E = -1.5 eV and E = -4 eV, (a)(b) (c)(d)(e) (f) Figure 2. The electronic band structures of TlSb without (a)and with (b) SOC. (c) 3D bulk Fermi surfaces and color-codedFermi velocities (red is high velocity). (d) Calculated Fermisurfaces cross section at the k z = 0 plane. (e) Calculatedenergy distribution at Dirac point in k x - k y plane. (a) (b)(d)(c) Figure 3. (a) Surface band structure for (001) plane alongprojected high symmetry points. The surface spectra of (001)plane with (b) E = E F , (c) E = -1.5 eV and (d) E = -4 eV. Figure 4. Wilson loops of six time-reversal invariant planesat (a) k = 0.0, (b) k = 0.5, (c) k = 0.0, (d) k = 0.5, (e) k = 0.0, (f) k = 0.5, where k , k , k are in units of thereciprocal lattice vectors. respectively. The surface band at E = -4 eV being deeplybelow E F can be ignored due to negligible contributionto the electronic properties of material. Due to TlSb hav-ing inversion-symmetry, its topology can be described byone strong topological index ν and three weak indices ν , ν , ν [41]. Thus we calculated the Wilson loops on sixtime-reversal invariant planes using WannierTools [39],and the results are shown in Fig. 4. According the thedefinition of Wilson loops [42, 43], the topological indicesare (1;000) indicating that TlSb is a strong topologicalmaterial.Second, we focus on the SC emerging in the Diracsemimetal TlSb discovered first by us. Figure 5(a) dis-plays the temperature dependence of resistivity, ρ ( T ),measured at zero field. With decreasing temperature, ρ xx decreases leisurely, exhibiting a poor metallic behav-ior, then drops to zero at 4.32 K, a superconductivitytransition occurring with the mid-temperature T midc =4.38 K, ∆ T c = 0.15 K. This superconducting transition isalso confirmed by the susceptibility measurement. Figure5(b) presents the temperature dependence of susceptibil-ity, χ ( T ), measured at H = 5 Oe with both zero-fieldcooling (ZFC) and field cooling (FC) process. It is clearthat a sharp diamagnetic transition emerges at 4.3 K,and the complete diamagnetism (4 πχ ∼ -1) below T c in-dicates the bulk superconductivity being from TlSb, sinceSb element is non-superconducting at ambient pressure.Figure 6(a) shows the field dependence of magnetiza-tion, M ( H ), measured at 2 K for a TlSb sample, exhibit- (cid:2) ( mW cm) T ( K ) (cid:1) ( mW cm) T ( K ) T c = 4 . 3 8 K p (cid:1) T ( K ) Z F C F C H = 5 O e T c = 4 . 3 K ( a ) Figure 5. (a) The temperature dependence of resistivity, ρ ( T ), of TlSb. The inset: enlarged ρ ( T ) near the super-conducting transition. (b) The temperature dependence ofmagnetic susceptibility, χ ( T ), measured at H = 5 Oe. ing a hysteresis, which indicates that TlSb is a type-IIsuperconductor. Then we measured the M ( H ) at var-ious temperatures below 4.5 K, as shown in Fig. 6(b).The lower critical field H c ( T ) can be estimated by thefield, at which M ( H ) curve starts to deviate from the lin-ear relationship. The obtained H c ( T ) is shown in theinset of Fig. 6(b), the lower critical field at zero temper-ature H c (0) = 148 Oe was obtained by the fitting usingthe GL relationship: H c ( T ) = H c (0)(1 − ( TT c ) ) (1)In order to get the upper critical field H c (0), wemeasured the superconducting transition temperature( T midc ) at various applied magnetic fields. As shown inthe inset of Fig. 7, the T c decreases, and the transitionwidth ∆ T c increases with increasing magnetic field. Byusing the GL formula H c ( T ) = H c (0)(1- t )/(1+ t ),where t is the normalized temperature t = T / T c , tofit the H c ( T ) data, the zero temperature upper criti-cal field H c (0) = 1.12 T was obtained, which is muchlower than the Pauli limit field H Pc (0) = 1.86 T c = 8.18 - 1 . 0 - 0 . 5 0 . 0 0 . 5 1 . 0- 3- 2- 10123 0 1 5 0 3 0 0- 1 . 6- 0 . 80 . 0 ( b ) M (emu/g) (cid:1) H ( T ) T = 2 K( a ) M (emu/g) H ( O e )
2 K H c1 (Oe) T ( K ) H c 1 ( 0 ) = 1 4 8 O e Figure 6. (a) Field dependence of magnetization M ( H ) mea-sured at 2 K. (b) The low field magnetization of TlSb atdifferent temperatures. The dashed line indicates the initiallinear magnetization curve. The inset shows the temperaturedependence of lower critical field, H c ( T ) determined by themagnetization curve deviating from linear. The red line is the H c ( T ) fitted by GL relation T. Then, the GL coherence length ξ GL (0) = 15.3 nmwas estimated by using the formula H c (0) = Φ /2 πξ GL ,where Φ is the quantum flux ( h /2 e ). The penetrationdepth λ GL (0) = 162 nm was estimated by using the for-mula H c (0) = (Φ /4 πλ GL (0))ln( λ GL (0)/ ξ GL (0)), andthe GL parameter κ GL = λ GL (0)/ ξ GL (0) = 10.6.We also measured the specific heat as a function oftemperature, C p ( T ), for TlSb in the temperature rangeof 0.5 - 5 K at both zero field and 3 T, respectively, asshown in Fig. 8. It is clear that the zero-field C p ( T ),compared with the C p ( T ) measured at 3 T ( > H c ,in this case bulk SC is completely suppressed), exhibitsa small and broad peak near T c , corresponding to thesuperconducting transition. No other anomaly was ob-served except for the peak near T c = 4.38 K, indicating T c m i d G - L F i t t i n g (cid:1) (T) T ( K ) (cid:1) ( mW cm) T ( K )
0 T1 . 5 T
Figure 7. The temperature dependence of upper critical field H c determined from resistivity measurements. The insetshows the low temperature ρ ( T ) curves measured at variousmagnetic fields. C p/ T (mJ mol-1 K-2) T ( K )
0 T 3 T D e b y e f i t C el/ g n T T / T c D C / g n T = 1 . 4 2 Figure 8. The temperature dependence of specific heat ofTlSb measured at 0 T (black circles) and 3 T (red circles)plotted as C p versus T . The solid blue line is a fit by theDebye model. The inset shows normalized electronic specificheat C el / γ n T versus T / T c at zero field. that no Tl impurities ( T c = 2.39 K) emerge in our sam-ple although a small amount of non-superconducting Sbimpurities was detected in the XRD. We fitted the lowtemperature C p ( T ) data measured at 3 T using the De-bye model: C p / T = γ n + β T + β T (2)where γ n is the Sommerfeld Coefficient, both the β T and β T are the phonon contributions to specific heat.The parameters γ n = 5.56 mJ mol − K − , β = 1.39mJ mol − K − , β = 0.44 mJ mol − K − were ob-tained. The inset of Fig. 8 shows the normalized ∆ C p γ n T TABLE I. Superconducting parameters of TlSbParameters unit value T c K 4.3 H c (0) Oe 148 H c (0) T 1.12 H Pc (0) T 8.18 ξ GL nm 15.3 λ GL nm 162 κ GL γ n mJ mol − K − C el / γ n T c = C p (0T) − C p (3T) γ n T as a function of the normalized tem-perature t = T / T c . The bulk superconducting temper-ature T c = 4.3 K was estimated by a entropy-balancemethod, consistent with the results from the resistivityand susceptibility measurements mentioned above. Thenormalized specific heat jump ∆ C el / γ n T c = 1.42 was es-timated, almost the same with the predicted value (1.43)by the Bardeen-Cooper-Schrieffer (BCS) theory [44], in-dicating that TlSb is a s -wave phonon-mediated super-conductor. The Debye temperature Θ D = 141 K wasestimated by using the formula Θ D = (12 π NR /5 β ) ,where N = 2 is the number of atoms in an unit celland the R = 8.314 J mol − K − is the molar gas con-stant. Using the obtained Θ D and T c values, we calcu-lated the electron-phonon coupling constant λ ep by usingthe McMillan formula [45]: λ ep = 1 .
04 + µ ∗ ln(Θ D / . T c )(1 − . µ ∗ ) ln(Θ D / . T c ) − .
04 (3)where µ ∗ = 0.13 is a typical value of the Coulomb re-pulsion pseudopotential for the intermetallic supercon-ductors. The obtained λ ep = 0.78 value is comparableto that of other superconductors such as PbTaSe ( λ ep = 0.74) [46] and Nb . Re . ( λ ep = 0.73) [47], suggest-ing TlSb is a moderately coupled superconductor. Theobtained superconducting parameters are summarised inTable I. IV. CONCLUSION
In summary, the calculations of the electronic bandstructure, the FS and the surface states show that TlSbwith a CsCl-type structure is a Dirac semimetal. Wemeasured the resistivity, magnetization, specific heat forthe polycrystalline TlSb sample. We first found thatTlSb is a type-II superconductor with T c = 4.38 K, H c (0) = 148 Oe, H c = 1.12 T and κ GL = 10.6. Thespecific heat results demonstrate it to be a moderatelycoupled s -wave superconductor. Although we can notdetermine which bands near E F to be responsible forSC, the coexistence of bulk SC with s -wave symmetryand the Dirac fermions on the surface in a single TlSbcompound provides an opportunity to realize Majoranazero energy mode. ACKNOWLEDGEMENTS
This research is supported by the National Key Pro-gram of China under Grant No. 2016YFA0300402and the National Natural Science Foundation of China(Grants No. NSFC-12074335 and 11974095) the Funda-mental Research Funds for the Central Universities, anopen program from the National Lab of Solid State Mi-crostructures of Nanjing University (Grant No. M32025). ∗ Corresponding author: [email protected][1] C. Beenakker and L. Kouwenhoven, Nat. Phys. , 618(2016).[2] M. Sato and Y. Ando, Rep. Prog. Phys. , 076501(2017).[3] Y. Maeno, H. Hashimoto, K. Yoshida, S. Nishizaki,T. Fujita, J. Bednorz, and F. Lichtenberg, Nature ,532 (1994).[4] G. M. Luke, Y. Fudamoto, K. Kojima, M. Larkin, J. Mer-rin, B. Nachumi, Y. Uemura, Y. Maeno, Z. Mao, Y. Mori, et al. , Nature , 558 (1998).[5] S. Kashiwaya, H. Kashiwaya, H. Kambara, T. Furuta,H. Yaguchi, Y. Tanaka, and Y. Maeno, Phys. Rev. Lett. , 077003 (2011).[6] A. P. Mackenzie and Y. Maeno, Rev. Mod. Phys. , 657(2003).[7] Y. S. Hor, A. J. Williams, J. G. Checkelsky, P. Roushan,J. Seo, Q. Xu, H. W. Zandbergen, A. Yazdani, N. P. Ong,and R. J. Cava, Phys. Rev. Lett. , 057001 (2010).[8] S. Sasaki, M. Kriener, K. Segawa, K. Yada, Y. Tanaka,M. Sato, and Y. Ando, Phys. Rev. Lett. , 217001(2011).[9] K. Matano, M. Kriener, K. Segawa, Y. Ando, and G.-q.Zheng, Nat. Phys. , 852 (2016).[10] S. Yonezawa, K. Tajiri, S. Nakata, Y. Nagai, Z. Wang,K. Segawa, Y. Ando, and Y. Maeno, Nat. Phys. , 123(2017).[11] Z. Liu, X. Yao, J. Shao, M. Zuo, L. Pi, S. Tan, C. Zhang,and Y. Zhang, J. Am. Chem. Soc. , 10512 (2015).[12] Y. Qiu, K. N. Sanders, J. Dai, J. E. Medvedeva, W. Wu,P. Ghaemi, T. Vojta, and Y. S. Hor, arXiv preprintarXiv:1512.03519 (2015).[13] S. Sasaki, Z. Ren, A. A. Taskin, K. Segawa, L. Fu, andY. Ando, Phys. Rev. Lett. , 217004 (2012).[14] L. Fu and C. L. Kane, Phys. Rev. Lett. , 096407(2008).[15] J. R. Williams, A. J. Bestwick, P. Gallagher, S. S. Hong,Y. Cui, A. S. Bleich, J. G. Analytis, I. R. Fisher, andD. Goldhaber-Gordon, Phys. Rev. Lett. , 056803(2012).[16] M.-X. Wang, C. Liu, J.-P. Xu, F. Yang, L. Miao, M.-Y.Yao, C. Gao, C. Shen, X. Ma, X. Chen, et al. , Science , 52 (2012).[17] F. Yang, F. Qu, J. Shen, Y. Ding, J. Chen, Z. Ji, G. Liu,J. Fan, C. Yang, L. Fu, and L. Lu, Phys. Rev. B ,134504 (2012).[18] E. Wang, H. Ding, A. V. Fedorov, W. Yao, Z. Li, Y.-F.Lv, K. Zhao, L.-G. Zhang, Z. Xu, J. Schneeloch, et al. ,Nat. Phys. , 621 (2013).[19] J. B. Oostinga, L. Maier, P. Sch¨uffelgen, D. Knott,C. Ames, C. Br¨une, G. Tkachov, H. Buhmann, and L. W.Molenkamp, Phys. Rev. X , 021007 (2013).[20] A. D. K. Finck, C. Kurter, Y. S. Hor, and D. J. Van Har-lingen, Phys. Rev. X , 041022 (2014).[21] M. Snelder, C. Molenaar, Y. Pan, D. Wu, Y. Huang,A. de Visser, A. Golubov, W. van der Wiel,H. Hilgenkamp, M. Golden, et al. , Supercond. Sci. Tech-nol. , 104001 (2014).[22] J. Wiedenmann, E. Bocquillon, R. S. Deacon,S. Hartinger, O. Herrmann, T. M. Klapwijk, L. Maier,C. Ames, C. Br¨une, C. Gould, et al. , Nat. Commun. ,10303 (2016).[23] D. Wang, L. Kong, P. Fan, H. Chen, S. Zhu, W. Liu,L. Cao, Y. Sun, S. Du, J. Schneeloch, et al. , Science ,333 (2018).[24] T. Machida, Y. Sun, S. Pyon, S. Takeda, Y. Kohsaka,T. Hanaguri, T. Sasagawa, and T. Tamegai, Nat. Mater., 1 (2019).[25] C. Chen, K. Jiang, Y. Zhang, C. Liu, Y. Liu, Z. Wang,and J. Wang, Nat. Phys. , 536 (2020).[26] X. Chen, M. Chen, W. Duan, H. Yang, and H.-H. Wen,Nano Lett. , 2965 (2020).[27] J. Yin, Z. Wu, J. Wang, Z. Ye, J. Gong, X. Hou, L. Shan,A. Li, X. Liang, X. Wu, et al. , Nat. Phys. , 543 (2015).[28] M. H. Fang, H. M. Pham, B. Qian, T. J. Liu, E. K.Vehstedt, Y. Liu, L. Spinu, and Z. Q. Mao, Phys. Rev.B , 224503 (2008).[29] L. Jin, X. Zhang, X. Dai, L. Wang, H. Liu, and G. Liu,IUCrJ , 688 (2019).[30] Y. Du, F. Tang, D. Wang, L. Sheng, E.-J. Kan, C.-G.Duan, S. Y. Savrasov, and X. Wan, npj Quantum Mater. , 1 (2017).[31] B. Predel and W. Schwermann, Z. Naturforsch. A. ,877 (1970).[32] A. C. Larson and R. B. Von Dreele, Report lAUR , 86(1994).[33] K. Schwarz, P. Blaha, and G. K. Madsen, Comput. Phys.Commun. , 71 (2002).[34] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996).[35] A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza,D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. , 2309 (2014).[36] N. Marzari and D. Vanderbilt, Phys. Rev. B , 12847(1997).[37] I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B , 035109 (2001).[38] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, andD. Vanderbilt, Rev. Mod. Phys. , 1419 (2012).[39] Q. Wu, S. Zhang, H.-F. Song, M. Troyer, and A. A.Soluyanov, Comput. Phys. Commun. , 405 (2018).[40] M. L. Sancho, J. L. Sancho, J. L. Sancho, and J. Rubio,J. Phys. F: Metal Phys. , 851 (1985).[41] L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. ,106803 (2007). [42] R. Yu, X. L. Qi, A. Bernevig, Z. Fang, and X. Dai, Phys.Rev. B , 075119 (2011).[43] A. A. Soluyanov and D. Vanderbilt, Phys. Rev. B ,035108 (2011).[44] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys.Rev. , 1175 (1957). [45] W. L. McMillan, Phys. Rev. , 331 (1968).[46] M. N. Ali, Q. D. Gibson, T. Klimczuk, and R. J. Cava,Phys. Rev. B , 020505 (2014).[47] A. B. Karki, Y. M. Xiong, N. Haldolaarachchige,S. Stadler, I. Vekhter, P. W. Adams, D. P. Young, W. A.Phelan, and J. Y. Chan, Phys. Rev. B83