Transition-Metal Pentatelluride ZrTe$_5$ and HfTe$_5$: a Paradigm for Large-gap Quantum Spin Hall Insulators

Quantum spin Hall (QSH) insulators, a new class of quantum matters, can support topologically protected helical edge modes inside bulk insulating gap, which can lead to dissipationless transport. A major obstacle to reach wide application of QSH is the lack of suitable QSH compounds, which should be easily fabricated and has large size of bulk gap. Here we predict that single layer ZrTe$_5$ and HfTe$_5$ are the most promising candidates to reach the large gap QSH insulators with bulk direct (indirect) band gap as large as 0.4 eV (0.1 eV), and robust against external strains. The 3D crystals of these two materials are good layered compounds with very weak inter-layer bonding and are located near the phase boundary between weak and strong topological insulators, which pave a new way to future experimental studies on both QSH effect and topological phase transitions.


I. INTRODUCTION
Topological insulators (TI) 1,2 can support gapless boundary states inside bulk insulating gap, which are topologically protected and robust against perturbations. Particularly in two-dimensional (2D) TIs, namely the QSH insulators 3,4 , the low energy (back) scattering of the edge states is prohibited by the time reversal symmetry, leading to the dissipationless transport edge channels and the QSH effect 5 . However, the known experimental verifications of QSH effect in HgTe/CdTe 6 and InAs/GaSb 7 quantum well structures require extreme conditions, such as the precisely controlled molecular-beam-epitaxy (MBE) and the ultra low temperature (due to the small bulk band gap of the order of meV), which greatly obstruct the further experimental studies and possible applications.
To be a "good" QSH insulator, the material must meet the following important criteria: (1) It must be a good layered materials to get the 2D system easily; (2) It must have large 2D bulk band gap to realize the QSH effect at high temperature. The first proposal for the QSH insulator in graphene 3 is practically useless due to its extremely small gap (10 −3 meV) induced by spin-orbit coupling (SOC) 8 , although graphene is so far the best 2D material which can be made easily even by scotch tape 9 . Recently Bi 2 TeI has been proposed to be another candidate for QSH insulator 10 , but its energy gap is also small. The Bismuth (111)bilayer is potentially a large-gap (about 0.2 eV) QSH insulator 11 , however, it is inter-layer strongly bonded and its fabrication is hard (not experimentally achieved yet). Several other proposals, such as the ultra thin tin films 12 and the Cd 3 As 2 quantum well 13 , have the same problem. Since they are not layered materials, the well controlled MBE technique is required to obtain the ultra thin film samples.
Here we report that simple binary ZrTe 5 and HfTe 5 , previously known as layered thermoelectric compounds, are inter-layer weakly bonded comparable to that of graphite. Their single layer sheets, which can be made in principle with no need for MBE, are QSH insulators with large bulk gap, and robust against the lattice distortions. Therefore, ZrTe 5 and HfTe 5 are the most promising 2D TI, which satisfy both the above conditions and pave a new way to more experimental studies on QSH effect. Moreover, our calculations show that their three dimensional (3D) crystals formed by the stacking of layers are located in the vicinity of a transition between strong and weak TI, which further makes it a perfect platform to study the topological quantum phase transitions.

II. COMPUTATIONAL METHODS
The ab initio calculations have been done by using the all electron full potential linearized augmented plane wave method implemented in the WIEN2k package 31 within the general gradient approximation (GGA). The choice of parameters, such as sampling of Brillouin Zone and cut-off of augmented plane waves, are carefully checked to ensure the convergence.
The spin-orbit coupling (SOC) is included self-consistently within the second variational method. The inter-layer binding energy, i.e., the unit area total-energy difference between the single layer sheet and the 3D bulk, is calculated by using the local density approximation (LDA) of Perdew and Wang type exchange correlation 32 . It is known that LDA can give the Van der Waals type inter-layer coupling energy reasonably well 33 . The maximally localized Wannier functions (MLWF) for Te p orbitals have been constructed by using home made code 34 based on OpenMX 35 , which can reproduce the band structures calculated from the first-principles quite well. The surface and edge states have been calculated for slabs by using a tight-binding model constructed from the MLWFs generated above.

A. Layered crystal structure
Both ZrTe 5 and HfTe 5 take the same structure and have very similar properties. They have been shown to possess interesting electrical transport properties, such as the giant resistivity anomaly 14,15 and the large thermopower 16,17 . They crystallize in the orthorhombic layered structure 18 with space group Cmcm (D 17 2h ) as shown in Fig. 1 (hereafter we take ZrTe 5 as an example). Trigonal prismatic chains of ZrTe 3 (marked by red dashed line) run along the a axis, and these prismatic chains are linked via parallel zigzag chains of Te atoms along the c axis to form 2D sheet of ZrTe 5 in the a-c plane. The sheets of ZrTe 5 stack along the b axis, forming layered structure. The prism of ZrTe 3 is formed by a dimer of Te d 1,2 atoms (with superscript d indicating the dimer and subscript 1, 2 numbering the atoms) and one apical Te a atom (with a meaning the apical), while the zigzag chain is formed by two Te z  ZrTe 5 shows strong quasi-2D anisotropy 19 . The prism chains and zigzag chains are connected through the apical Te atoms, and the Te-Te bond length between two chains is just about 0.4Å longer than that in the zigzag chain. If the Zr and Te d dimer atoms are neglected, the remaining apical Te a and zigzag Te z atoms can be viewed as a waved grid of Te square lattice sheet (see Fig.1(d)), leading to a stable quasi-2D structure. Each ZrTe 5 layer is nominally charge neutral, and the inter-layer distance (along b axis) is quite large (about 7.25Å), suggesting the weak inter-layer coupling of presumably Van der Waals type.
We have calculated the inter-layer binding energies for different compounds by using the first-principles total energy method (see method). The results shown in Fig. 2 suggest that the inter-layer binding energy of ZrTe 5 or HfTe 5 is as weak as that of graphite, and is much smaller than that of Bi 2 Se 3 and Bi (111) graphene from graphite simply by scotch tape 9 , the comparably weak inter-layer binding energy suggests that the single layer of ZrTe 5 (or HfTe 5 ) may be formed in a similar simple and efficient method.

B. 2D Quantum Spin Hall Insulator
We now focus on the single layer of ZrTe 5 in the a-c plane. Its structure is fully optimized by theoretical calculation. It is found that the relaxed lattice constants (a=4.036Å and c=13.843Å) and the internal atomic coordinates are all very close to their 3D bulk experimental values (less than one percent difference) 18 , again suggesting the stability of the single layer sheet. The calculated density of states shown in Fig. 3 For the convenience of discussions, hereafter we define our coordinate system with x, y being along the a and c-axis, respectively, and z being out of the 2D plane. We choose the origin of the coordinate system to be located on the Zr site. The calculated band structures for a single layer ZrTe 5 are shown in Fig. 3 (b) and (c). In the case without SOC, the system is a semimetal with a band-crossing along the Γ-X direction, implying the existence of band inversion around the Γ. The band crossing is unavoidable because the two bands belong to different representations distinguished by the m xz symmetry along the  20 , and conclude that it is a non-trivial QSH insulator with Z 2 =1. Since the band gap is mostly determined by SOC strength, we wouldn't expect much error bar of GGA type calculations (as it usually does for conventional semiconductors). Indeed, we have checked the electronic structure by using the hybrid functional HSE06, 21 as well as the modified topology unchanged.
The 2D non-trivial insulating state in single layer ZrTe 5 should support topologically protected conducting edge states. Two idea edges along the x-and y-direction are studied by tight-binding slab models (see method). The termination is between the prism chain and the zigzag chain for the former, and is between the two Te z atoms of the zigzag chain for the later. The results are shown in Fig. 3(d) and (e). For the case of cutting along the x direction, the left and right edges of our slab are not symmetric, which leads to two separated Dirac cones atΓ: one for the edge with prism chain termination (red), and other one for the zigzag chain termination (blue). While for the case of cutting along the y direction, both sides of the slab are symmetric and two Dirac cones (located at opposite sides) are energetically degenerate. There is another difference between x-and y-edge. The Dirac cones are atΓ in the former, while atȲ in the later.

C. Band inversion mechanism
The band inversion in single layer ZrTe5 is not due to SOC, instead it is mainly due to the non-symmorphic features of the space group. Therefore the physics of band inversion here can be understood without SOC, and the only effect of SOC is to open a energy gap afterwards. We note that the space group of the system is P mmn (D 13 2h ), which is a nonsymmorphic one with the inversion center located not at the origin but at (1/4,1/4). As an important consequence of such space group symmetry, the eigenstates at all zone boundary TRIM points (namely, X, Y, and M) are all four-fold degenerate (including spin degrees of freedom) with two of them having even parity and the other two odd parity. This is because that at all the three zone boundary points the inversion operation anti-commutes with two important mirror operators m xz and m yz . (see detail in Appendix). The equal number of parities guarantees that any band inversion at those zone boundary TRIM points will not change the topology of the system and the Z 2 index of the material is fully determined by the energy order of the bands at the Γ point, where only the Kramer degeneracy holds. We further notice that the mirror symmetry m yz will transfer all the atoms into themselves. Therefore, the atomic p orbitals can be classified into two classes, namely, the p x orbital and the p y/z orbital, because they will not mix in the eigenstates of Γ point (in the absence of SOC).

Our calculations suggest that the band inversion at Γ happens between the zigzag chain
Te z -p x and the prism chain dimer Te d -p y states, as shown schematically in Fig. 4. At the atomic limit, they are all four-fold degenerated since there are four equivalent Te atoms for each class. The strong intra-chain covalence bonding will split them into double degenerate bonding and anti-bonding states, respectively, and then the weak inter-chain coupling will further split them into singly degenerate states. It is important to note the symmetry difference between the two chains. The two dimer Te d 1,2 atoms in the same prism chain are related by the mirror m xz symmetry, while the two Te z 1,2 atoms in the same zigzag chain are related by the inversion symmetry. The bonding and antibonding states caused by the intrachain coupling can therefore be distinguished by the eigenvalues m xz =±1 for the former, and the parity p=±1 for the later. Due to the strong intra-chain coupling, the band inversion is introduced between the bonding Te d states with m xz =1 and the antibonding Te z states with p=-1, as illustrated in Fig. 4. Within the two bonding Te d states of m xz = 1, only one state has odd parity. As a result, its occupation changes the total parity of the occupied states at Γ, which leads to the QSH state.

D. Strain effects
We emphasize that, unlike the situation in typical TI such as the Bi 2 Se 3 family compounds 23 , the band inversion in ZrTe 5 single layer is not due to SOC, instead it is due to the special non-symmorphic space group features as discussed above. As long as the inter-chain coupling is not strong enough to reverse the band ordering again (between the (m xz =-1, p=+1) state and the (m xz =1, p=-1) state, as shown in Fig. 4), the QSH state should be stable. We have checked the stability of this QSH state as function of strains. Without lost of generality, here we assume that the lattice parameters a and c are equally scaled by possible external potentials, and for each fixed volume we fully optimize the internal coordinates. The results shown in Fig. 5 suggest that the QSH state survives even if the volume is expanded by more than 20% or compressed by 10%, indicating its robust stability against the strains. This makes it highly adaptable in various application environments. In the case of compressing, the too much enhanced inter-chain coupling will finally kill the QSH state, as discussed above.  Energy(eV)

intra-dimer
Energy(eV) (see main text for details).

E. 3D weak and strong TI
Given the QSH state in single layer ZrTe 5 , it is important to check if the stacked 3D compound is a weak or a strong TI 24 . Interestingly, if we use the experimental lattice parameters 18 , we get a 3D strong TI with topological invariant (1; 001); however, if the optimized lattice parameters, which are only about 2% larger than experimental ones, are used, we get a 3D weak TI state 25 with topological invariant (0; 001). The difference between two solutions are due to the tiny overlap of the states at the Γ point of 3D BZ, which inverse the parity (see Fig. 6). If we use the mBJ potential 22 instead of GGA for the above two calculations, we get both strong TI solutions. We have to conclude that the required accuracy to distinguish these two states for reality is beyond the ability of our present calculations. Nevertheless, our results suggest that the 3D compound is located very close to the boundary between the weak and strong TI, which provides us an excellent opportunity to study the topological phase transition experimentally. In Fig. 6 The breaking of time reversal symmetry (TRS) will allow us to study the possible quantum anomalous Hall effect 26,27 . Moreover, due to the strong 2D features of these materials, the edge states will appear not only at the edge of the single-layer samples but also at the terrace edges on the top surface of the 3D materials. When such terrace edges are covered by a s-wave superconducting layer, 1D topological superconductors with Majorana bound states at the end of the terrace edges will be induced by the superconducting proximity effect [28][29][30] .

V. ACKNOWLEDGEMENTS
We acknowledge the supports from the NSF of china and the 973 program of China (No. 2011CBA00108 and 2013CBP21700).

VI. APPENDIX
In this appendix, we prove that for single layer ZrTe 5 the energy bands at all the zone boundary TRIM points are at least doubly degenerate with opposite parity. For simplicity, we neglect the spin degree of freedom and the conclusion can be easily generalized to the spin-1/2 case. Since the energy bands of ZrTe 5 near the Fermi level are formed by Te 5porbitals, the Wannier functions centered at each atomic sites with particular atomic feature (like p x , p y and p z ) carry a natural representation for the space group. First we define the following Wannier basis in the momentum space, , where i is the unit cell index, µ denotes the different atomic positions in a unit cell and α labels different orbitals on a same atomic position. When a general space group operator {P R | t R }, with P R and t R representing the rotation and translation part of the operator respectively, acts on the Wannier basis, we have where we have take t R + Rτ µ = R µ 0 + τ µ and α = Rα, Z R and O R are matrices describing the transformation in atomic positions and orbitals respectively. We would emphasize that in the above equation φ µ αRk (r) does not necessary equal to φ µ αk (r), even when Rk = k + G n with G n representing a vector in reciprocal lattice. They can differ by a phase factor due to the specific choice of the gauge fixing condition here, which contains a phase factor e ikτµ for orbitals located at different atomic position τ µ . Using the Wannier representation introduced above and choosing the Zr site as the origin of the coordinate frame, we pick two operators, inversion and mirror m yz , which can be written as the following between the above two is the phase factor, which differs by −1 for (π, 0) and (π, π). Therefore for these two points we have {I| 1 2 , 1 2 } · {m yz |0, 0} = −{m yz |0, 0} · {I| 1 2 , 1 2 }. Because both the above symmetry operators commute with the Hamiltonian, we can easily prove that if φ k is eigen state of Hamiltonian with given parity, the state {m yz |0, 0}φ k is another eigen state of Hamiltonian but with opposite parity. Therefore we have proved that at (π, 0) and (π, π), all the states come in pairs with opposite parity. Replacing m yz by another mirror m xz , we can easily prove that another zone boundary TRIM point (0, π) has the same property.
In conclusion, at all these zone boundary TRIM points, the eigen states always come in degenerate pairs with opposite parity.

VII. AUTHOR CONTRIBUTIONS
H.M.W. did the material search and the first principle calculations. X.D. and Z.F. did the symmetry and topology analysis. All the authors contributed to the analysis of the computational data and writing the manuscript. * Electronic address: daix@aphy.iphy.ac.cn