Topological insulators in the NaCaBi family with large SOC gaps

By means of first-principles calculations and crystal structure searching techniques, we predict that a new NaCaBi family crystallized into the ZrBeSi-type structure (\ie $P6_{3}/mmc$) are strong topological insulators (STIs). Taking $P6_{3}/mmc$ NaCaBi as an example, the calculated band structure indicates that there is a band inversion between two opposite-parity bands at the $\Gamma$ point. In contrast to the well-known Bi$_2$Se$_3$ family, the band inversion in the NaCaBi family has already occured even without spin-orbit coupling (SOC), giving rise to a nodal ring surrounding $\Gamma$ in the $k_z=0$ plane (protected by $M_z$ symmetry). With time reversal symmetry $(\cal T)$ and inversion symmetry $(\cal I)$, the spinless nodal-line metallic phase protected by $[{\cal TI}]^2=1$ is the weak-SOC limit of the spinful topological insulating phase. Upon including SOC, the nodal ring is gapped, driving the system into a STI. Besides inversion symmetry, the nontrivial topology of NaCaBi can also be indicated by $\bar{6}$ symmetry. More surprisingly, the SOC-induced band gap in NaCaBi is about 0.34 eV, which is larger than the energy scale of room temperature. Four other compounds (KBaBi, KSrBi, RbBaBi and RbSrBi) in the family are stable at ambient pressure, both in thermodynamics and lattice dynamics, even though the gaps of them are smaller than that of NaCaBi. Thus, they provide good platforms to study topological states both in theory and experiments.


I. INTRODUCTION
Topological insulators (TIs) have been studied broadly and deeply in both theory and experiments, since the two-dimensional quantum spin Hall effect (2D TI) was proposed in 2005 1 . By generalizing the concept of TIs, the three-dimensional (3D) TIs were predicted in the Bi 2 Se 3 family 2 , which were verified in the experiments later [3][4][5] . TIs embody a new quantum state characterized by a topological invariant rather than a spontaneously broken symmetry. The nontrivial topology of TIs is well defined by the wave functions of occupied bands, as long as there is a continuous energy gap between the conduction bands and valence bands in the whole Brillouin zone (BZ). Then, the T -invariant TIs are generalized to topological crystalline insulators (TCIs) [6][7][8][9][10] , in which the lattice symmetries play important roles. Generally, symmetry eigenvalues and irreducible representations (irreps) can indicate topological property. For example, the wellknown Fu-Kane parity indices are widely used for centrosymmetric TIs 11 , and the recently established symmetry indicators and topological quantum chemistry are very efficient to diagnose the nontrivial topology of TIs and TCIs [12][13][14][15][16][17][18][19][20] . Although many candidates for TIs/T-CIs have been proposed, the TIs protected by sixfold rotoinversion symmetry [i.e., the sixfold rotoinversion 6 (≡ IC 6 ) is a sixfold rotation (C 6 ) followed by inversion (I)] are rarely reported.
The 3D TIs have attracted tremendous attention because their surface states possess odd spin-momentum locked Dirac cones, thus, the backscattering channels of the surface states are suppressed, leading to dissipationless electronic transport 2,3,21,22 . In proximity to an swave superconductor (SC), the topological surface states can acquire an effective p-wave superconducting paring, resulting in a 2D topological SC 23 . In addition, the TIs can be tuned into many interesting topological states. For instance, the quantum anomalous Hall (QAH) insulator is found in the Cr-doped Bi 2 Te 3 24,25 . Furthermore, the compound MnBi 2 Te 4 is predicted to be an antiferromagnetic TI at zero magnetic field and a Weyl semimetal under the out-of-plane magnetic field 26 . Recently, the QAH effect has been observed in few-layer thin films 27 . Moreover, the unconventional SC has been proposed and some experimental evidences were revealed in Cu-doped Bi 2 Se 3 28,29 . The strong evidences of Majorana zero-energy mode within a vortex in the TI/SC heterostructure -Bi 2 Te 3 /NbSe 2 -have been proposed recently 30 . Therefore, the TI candidates are of great interest, especially those with relatively large band gaps, which not only benefit the the future applications of spintronics and quantum computation from the dissipationless surface states, but also provide good platforms to study the interplay between the topological phase and symmetry-breaking phases. Though many TI candidates and corresponding interesting properties have been proposed 2,21,22,31,32 , TIs with relatively large band gaps and clear 2D Dirac-cone surface states as the Bi 2 Se 3 family are rare.
In this work, using crystal structure searching techniques and first-principles calculations, we predict that  1. (a) Calculated enthalpy relative to that of P 63/mmc phase vs. pressure. (b) Phonon dispersions for the P 63/mmc phase at the ambient pressure, which justify the lattice dynamical stability of this phase. (c) Five crystal structures crystallized into the P 63/mmc, P 3m1, P nma, F43m and P42m phases of stoichiometric NaCaBi near the convex hull from the formation enthalpies.
NaCaBi can be stable in the space group (SG) of P 6 3 /mmc at pressures up to 10 GPa from both the thermal dynamics and the lattice dynamics. Calculations of the band structures indicate that without SOC, P 6 3 /mmc NaCaBi is a topological semimetal with a nodal ring surrounding the Γ point in the k z = 0 plane. Upon including SOC, the crossing points of the nodal line are all gapped and the system becomes a STI. We find that the nontrivial topology can be indicated not only by inversion symmetry, but also by the6 symmetry. Moreover, it is exciting that the topologically nontrivial band gap is about 0.34 eV due to strong SOC, which is even larger than the energy scale of room temperature. At last, four other dynamically stable XYBi compounds with similar electronic structures and STI natures are proposed in Section 4 of Appendix.

A. The crystal structures of stoichiometric NaCaBi
We used crystal structure prediction techniques integrated in USPEX [33][34][35] and the machine learning accelerated crystal structure search method 36 to find the best candidates of NaCaBi under pressure. From the phase diagrams based on convex hull analysis of formation enthalpies for the Na-Ca-Bi ternary system at 0 GPa and 10 GPa shown in Fig. 5, we find stoichiometric NaCaBi crystallizes in P nma phase at 0 GPa and P 6 3 /mmc phase at 10 GPa, respectively. We extract several structures of stoichiometric NaCaBi near the two stable phases from the formation enthalpy, which are in SGs of P 6 3 /mmc, P 3m1, P nma, F43m and P42m, respectively. The enthalpy-pressure ( H − P ) curves between these phases plotted in Fig. 1(a) exhibit the best ones from our structural predictions at pressures up to 10 GPa. At ambient pressure, the most stable structure of NaCaBi is found to be P nma phase, and the enthalpy of P 6 3 /mmc phase is just 0.17 eV/f.u. higher than that of P nma phase. Our calculations reveal that P 6 3 /mmc phase becomes preferred at pressures higher than 1 GPa.
As shown in Fig. 1(b) and Fig. 8, there is no phonon mode with negative frequency in the phonon spectra of P 6 3 /mmc NaCaBi under pressures from 0 GPa to 10 GPa, which indicates that the metastable P 6 3 /mmc phase could be quenched to ambient pressure if they can be synthesized at higher pressure. These results are in accordance with the earlier work 37,38 . As shown in Fig. 7(c), calculations of the ab initio molecular dynamic (AIMD) simulations at ambient pressure and T = 600 K indicates that there is no structural collapse after 10 ps (10000 steps). In addition, motivated by earlier works 39,40 , how the parameters evolves under pressure, strain potential and the possible substrates are also discussed in Section 8 of Appendix. As shown in Fig. 8 and Table VI, we also list many other XYBi compounds in the P 6 3 /mmc phase which are dynamically stable at low pressures. Besides NaCaBi, other four XYBi compounds (i.e., KBaBi, KSrBi, RbBaBi, and RbSrBi) have similar band structures, thus the same STI nature with NaCaBi (see Section 4 of the Appendix). Different from NaCaBi, the P 6 3 /mmc phase of KBaBi, KSrBi, RbBaBi and Rb-SrBi are always energetically preferred from 0 GPa to 10 GPa, as shown in Fig. 6. Moreover, the NaCaBi system possesses the largest inverted band gap, as shown in Table I. In the following, we will focus on NaCaBi with P 6 3 /mmc phase in the main text. From the band structures of P 6 3 /mmc NaCaBi without SOC shown in Fig. 2(a), we can find two band crossings along K-Γ and Γ-M at E F . The orbital-resolved band structures in Fig. 2(a) indicate that the valence band maximum is contributed by the Bi-p states, while the conduction band minimum is from the Ca-s states (which also hybridize with the Bi-s, Na-s and Ca-d states). Calculations of the irreps 41 indicate that the two bands belong to GM1+ and GM2− irreps respectively [ Fig. 2(a)], in the convention of the Bilbao Crystallographic Server notations at Γ (GM). The parities of the occupied bands at eight time-reversal invariant momenta (TRIMs) are given in Table II. Thus, we conclude that the band inversion appears between opposite-parity bands around Γ without SOC, which guarantees the existence of nodal line(s) in the spinless systems with T and I. Furthermore, this nodal ring surrounding Γ lies in the k z = 0 plane due to the presence of M z symmetry (i.e., the two inverted bands have different M z eigenvalues).
The spinful topological insulating phase can be considered as a (strong-)SOC limit of the spinless nodal-line metallic phase protected by [T I] 2 = 1 12 . As shown in Fig. 2(c), we can find that SOC induces a visible band gap between the two inverted bands around Γ, but it does not change the ordering of energy bands at Γ. The nodal ring is fully gapped, thus, the NaCaBi becomes a STI. More surprisingly, the SOC-induced band gap is as large as 0.34 eV, which is comparable with thermal energy of room temperature. The P 6 3 /mmc phase of NaCaBi family compounds provides good platforms to study topological states both in theory and experiments. Considering the possible underestimation of the band gap by Perdew-Burke-Ernzerhof (PBE) method 42 , the band inversion can be further confirmed by the calculations using hybrid Heyd-Scuseria-Ernzerhof (HSE) functional 43 . The results shown in Section 5 of the Appendix indicate that the SOC-induced band gap is almost unchanged, although the inverted gap at Γ (i.e., Symmetry eigenvalues (or irreps) play key roles in identifying TIs/TCIs in solids [14][15][16] . As we know, SG 194 belongs to a Z 12 classification. Using irvsp 41 and Check Topological Mat 20 , all symmetry indicators are calculated to be z 2w,i=1,2,3 = 0; z 4 = 3; z 6m,0 = 5; z 12 = 11. By performing the symmetry analysis, we find that there are two essential symmetries indicating the nontrivial topology in this system. One is the inversion symmetry, while the other is the6 symmetry. In other words, even though the inversion is broken when the SG is reduced to SG 174 (P6), the nontrivial topology can be still indicated by z 3m,0 = 1 and z 3m,π = 0. With M z ≡ [6] 3 , these indicate that the mirror Chern numbers of k z = 0 and k z = π planes are 1 and 0, respectively. Thus, this system has to be a 3D STI.
C. The surface states of P 63/mmc NaCaBi Exotic topological surface states serve as significant fingerprints to identify various topological phases. Based on the tight-binding model constructed with the maximally localized Wannier functions (MLWFs) and surface Green's function methods [44][45][46] , we have calculated the corresponding surface states of this system with SOC to identify its topological properties (see Section 6 of the Appendix). As shown in Figs. 3(a) and 3(b), the 2D Dirac cone is obtained on (010) and (001) surfaces. Thus, the calculated surface states agree well with the analysis of the symmetry indicators, which confirm the STI nature of NaCaBi.

D.
Low-energy effective model We find that the topologically nontrivial nature is determined by the low-energy dispersions near E F around Γ, which motivates us to build the effective k·p Hamiltonian to get further insights in this system. In the presence of SOC, the irreps for the lower energy bands at Γ are Γ + y . This is nothing but the well-known 3D TI Hamiltonian in the Bi 2 Se 3 family 2 in the case of M ·B 1,2 > 0, the band inversion happens between the two different-parity bands. By fitting the energy dispersions, we obtain the parameters: As aforementioned, the nontrivial topology of this system can be also protected/indicated by6 symmetry. If we introduce the additional δ 1 (δ 2 ) term as shown in Eq. (2), C 2x (C 2y ) and M y (M x ) will be broken. Therefore, by simply adding the two additional terms (both of which break C 6z and I), the6-invariant Hamiltonian is obtained below, In the presence of6 symmetry, the C 3z and M z symmetries are still preserved. Thus, one can check that in the case of M · B 1,2 > 0, the mirror Chern number of the k z = 0 plane is 1. On the other hand, one can also compute the symmetry indicators of SG 174 (P6), which is founed to be z 3m,0 = 1. This implies that the Z 2 invariant for the k z = 0 plane is 1. Combined with the trivial Z 2 invariant for the k z = π plane, one can conclude that z 3m,0 = 1 and z 3m,π = 0 yield the STI nature of P 6 3 /mmc NaCaBi. Suppose that we move the two Na atoms at 2a Wyckoff positions reversely along z axis a little (e.g., two percent of the third lattice vector along z axis), this distortion preserves6 symmetry, but breaks inversion symmetry and C 6z symmetry. In this disturbed structure, the δ 1 and δ 2 are fitted to be 0.4 eVÅ 2 and 0 eVÅ 2 . The STI phase is indicated by the symmetry indicators: z 3m,0 = 1 and z 3m,π = 0, since there is no band crossing occurring in this process.

III. DISCUSSION
In carrier-doped TIs, odd-parity pairing and nematic pairing have been theoretically proposed to realize topological superconductivity 48,49 . Some supporting evidences in Cu-doped Bi 2 Se 3 are revealed in experiments but the pairing symmetry is still under active debate [50][51][52][53][54] . In analogous to Bi 2 Se 3 , superconductivity and topological paring states may be achieved with electron doping in NaCaBi. One prominent feature of P 6 3 /mmc NaCaBi is that the band inversion is kept with SOC and is significantly large, resulting a double-peak valence band around Γ, as shown in Fig. 2(c). Within a relatively large hole doping region, the Fermi surface is a torus around the Γ point rather than spheres in other TIs, as shown in Fig. 2(d). As the topology of the Fermi surface is distinct from electron-doped case, the pairing state can be also different and may be unconventional. For example, for the intra-orbital spin singlet pairing c +↑ c +↓ − c −↑ c −↓ , which belongs to a trivial irrep, two nodal lines of gap functions in band space are expected on the torus. We leave the question of determining possible pairing states by solving effective models to future work. Nevertheless, we investigate the superconducting properties by electron-phonon coupling calculations and the Allen-Dynes modified McMillian equation 55 . Considering 0.06 hole dopping in NaCaBi system at ambient pressure, Tc is estimated to be about 1.5 K with a commonly used screened Coulomb potential µ * = 0.10. Most of the contributions to electron-phonon coupling constant arise from the B 2g vibrational mode, i.e., vibrations of Ca atoms and Bi atoms along z axis, as shown in Fig. 4(b).
In summary, we propose a new NaCaBi family which are STIs with large SOC-induced band gaps, among which KBaBi, KSrBi, RbBaBi and RbSrBi are stable at ambient pressure, both in thermodynamics and lattice dynamics. In terms of the NaCaBi system, the corresponding symmetry analysis indicates that there exists a clear band inversion between two opposite-parity bands near E F . Without SOC, the band inversion gives rise to

Ab initio calculations.
We performed firstprinciples calculations based on the density functional theory (DFT) using projector augmented wave (PAW) method implemented in the Vienna ab initio simulation package (VASP) 56 for the calculations of the structure optimization and electronic structures. The generalized gradient approximation (GGA), as implemented in the Perdew-Burke-Ernzerhof (PBE) functional 42 was adopted. The cutoff parameter for the wave functions was 850 eV. The BZ was sampled by Monkhorst-Pack method 57 with a k-spacing of 0.025 × 2πÅ −1 for the 3D periodic boundary conditions. The phonon calculations were performed by finite displacement method implemented in the PHONOPY code 58 , with a 2 × 2 × 2 supercell. A kmesh of 20×20×16 was used to calculate the Fermi surface of the valence band near the Fermi level. Electron-phonon coupling calculations are performed in the framework of Density functional perturbation theory, as implemented in the quantum-espresso code 59 .

Crystal structure searching
We perform a variable composition calculation to get the phase diagrams of Na-Ca-Bi ternary system at ambient condition (0 GPa) and 10 GPa. At ambient condition, we find P nma NaCaBi is the unique ternary phase lying in the convex hull (denoted by the solid red circle in the center of the equilateral triangle), as shown in Fig. 5(a). Meanwhile, there exist three other metastable ternary candidates, i.e., P4/mmm-NaCaBi 2 , Cmmm NaCa 2 Bi 3 and R3m Na 4 CaBi 2 , above the convex hull (denoted by the hollow red circles in the the equilateral triangle). The corresponding enthalpies of these metastable ternary candidates above convex hull (denoted by H) of P 4/mmm 1. The "SG" column denotes the space group in which NaCaBi crystallises into. 2. {α, β, γ} and {a,b,c} denote the lattice parameters of the conventional cell (while not primitive cell) in each phase. 3. The "Wyckoff positions" column gives the Wyckoff positions of Na, Ca and Bi in each phase. "1a ⊕ 1c", "1b ⊕ 1c" and "1a ⊕ 1b" in P 3m1 phase indicates that there exist two nonequivalent "'Na"/"'Ca"/"'Bi" atoms occupying (1a and 1c)/(1b and 1c)/(1a and 1b) Wyckoff positions in this phase, respectively. NaCaBi 2 , Cmmm NaCa 2 Bi 3 and R3m Na 4 CaBi 2 are 16.0 meV/atom, 13.3 meV/atom and 19.7 meV/atom, respectively. We have listed the structure parameters of these ternary structures in Table. IV.
At 10 GPa, we find P 6 3 /mmc NaCaBi and Cmmm NaCa 2 Bi 3 (denoted by two solid red circles lying inside the equilateral triangle) are two stable ternary structures lying in the convex hull, as shown in Fig. 5(b). In addition, the convex hull analysis indicates that there are three other structures (C2/m Ca 3 Bi 2 , Cmmm NaCaBi 2 and P 4/mmm NaCa 2 Bi) lying above the convex hull, which indicates they are metastable candidates. The corresponding structure parameters are listed in Table. V. It should be noted that C2/m Ca 3 Bi 2 and Cmmm NaCaBi 2 are metastable candidates with only 0.64 meV/atom and 1.2 meV/atom above the convex hull. We also remind that binary Na-Bi and binary Ca-Bi phases in the convex hull are almost in accordance with earlier works 60,61 .
We do not perform ternary structural searching for the other four XYBi candidates, because it will be time consuming. But, similar with NaCaBi, the thermal stability between the P 6 3 /mmc phase and the P nma phase are discussed. As shown in Fig. 6, the P nma phase of NaSrBi is always preferred under pressures from 0 GPa to 10 GPa, while for cases of KBaBi, KSrBi, Rb-BaBi and RbSrBi, the P 6 3 /mmc phase are always preferred. It means the STI phase can be obtained in the KBaBi, KSrBi, RbBaBi and RbSrBi systems without any quenched progress.
collapse was observed after 10 ps (10000 steps), which indicates the thermal stability of NaCaBi at ambient pressure.  Structures with "*" behind denote those possessing similar band structures with P 6 3 /mmc NaCaBi.
RbMgBi) and the STI phase. Since we are focused on the STI phase, several of these dynamically stable compounds of the STI phase, having the similar band structures with P 6 3 /mmc NaCaBi in the main text, are denoted by "*" symbols.

HSE calculations
It is known that the PBE functional tends to underestimate the band gap, we recalculated the band structures of P 6 3 /mmc NaCaBi system with the HSE06 functional 43 , as shown in Fig. 9. We find the energy dispersions for low energy bands using the HSE06 functional is similar to that using PBE functional, which indicates that the topological nature is solid.

Wannier-based tight-binging model and surface states
Surface state plays an important role on justifying/characterizing TIs from the bulk-boundary correspondence. Here we choose Ca-s, Ca-d and Bi-p as the projected orbits to build the Wannier-based tightbinging model, it reproduces the DFT band structures well, as shown in Fig. 10. Then, we calculated the surface states based on this Wannier-based tight-binging model of NaCaBi on the (010) and (001) surfaces, as presented in Fig. 3 in the main text.

Matrix presentations of symmetry operators
Under the basis of Γ + 7 and Γ − 7 , the symmetry operators (e.g., T , I, C 6z , C 2x , C 2y ) are represented below: where K denotes the complex conjugation, τ x,y,z (σ x,y,z ) are Pauli matrices in the orbital (spin) space. Their (O = P, Q, R, T ) representations in momentum space are given The low-energy effective k · p Hamiltonian satisfies the condition placed by the symmetry O: After considering all symmetry restrictions of D 6h , the low-energy effective model is derived as given in the main text. From the electronic properties of P 6 3 /mmc NaCaBi at 0 GPa, 5 GPa and 10 GPa shown in Fig. 8, we find that there is an electron pocket of the conductance band around K dropping to the Fermi level quickly under pressure. And when the applied pressure is as large as 10 GPa, this pocket goes across the Fermi level, which leads P 6 3 /mmc NaCaBi to be a metal. As shown in Fig. 11(a), we can find the lattice parameters a(= b) and c decrease under pressure. Moreover, the lattice parameters a/c decrease about 1%/1.3% at 1 GPa, from which the P 6 3 /mmc phase of NaCaBi becomes preferred. It should also be noted that though the P 6 3 /mmc NaCaBi changes from an insulator to a metal from 0 GPa to 10 GPa, the inverted direct gap is always kept. Thus, the STI nature with topological surface states is kept up to 10 GPa.
Looking along the c axis, the P 6 3 /mmc NaCaBi is a layered hexagonal structure with its layer lying in the xy plane. Thus, we can introduce substrates perpendicular to the layer plane. As shown in Table VII, we have introduced 31 substrates with the lattice mismatch less than 1% along the x and y axes to assist growth teams attempting to synthesize the P 6 3 /mmc phase of NaCaBi. Some commonly used crystals such as KTiO 3 and BaNiO 3 only have lattice mismatch with NaCaBi of about -0.1% and 0.3%, which are very small and can be good substrate candidates. However, it should be noted that, as theoreticians, our proposal is very rough and more complicated experimental details should be considered before real experiments.
Similar with applying pressure, imposing various strains are also considered as convenient methods to tune the electronic properties in the condensed systems. Here we focus on applying biaxial linear strain along a and b axes which keeps the symmetries of this system unchanged. From Fig. 11(b), we find the strain energy of this kind of biaxial linear strain increase parabolically with −5% ≤ ε ≤ 5%, which indicates this kind of strain can be seen as elastic strain in a very large range. Thus, it offers the experimental team large degree of freedom to choose potential substrates candidates. 2. > 0 denotes the lattice parameter a of the substrate candidate is larger than that of P 6 3 /mmc NaCaBi, while < 0 denotes the lattice parameter a of the substrate candidate is smaller than that of P 6 3 /mmc NaCaBi.