Prediction of wide-gap topological insulating phase in metastable BiTeI

BiTeI crystals possess many unique pressure-induced electronic properties; however, high-pressure conditions limit their practical applications to some extent. We conducted crystal structure prediction from 3000 structures to determine several structures with high thermodynamic stability. By studying the topological invariants, we found one strong topological insulator with inversion symmetry and one weak topological insulator. Spin splitting was found in two structures, and their spin texture was analyzed by effective SOC Hamiltonian and symmetry. Our study showed that some structures are not only topologically nontrivial but can also simultaneously possess Rashba-type spin splitting with a band gap of up to 860 meV.

. Here, L and S are the orbital and spin angular momentum operators respectively. Z is the atomic number, c is the speed of light, n and l are the principal and orbital quantum numbers respectively. Thus, the strength of atomic SOC depends on the fourth power of the atomic number; therefore many of the current spintronic materials are found in systems with heavy elements, such as Bi, 8) Sb, 9) Pb, 10) etc. Bismuth compounds have attracted attention for spintronics applications since Bi 2 Te 3 was experimentally confirmed as the first topological insulator. 11) Recently, layered bismuth tellurium iodide (BiTeI) was found to possess considerable Rashba-type spin splitting, 12,13) making it a potential material for spintronics. M.S. Bahramy et al. predicted that when BiTeI was under pressure, band inversion and topological phase transitions were induced owing to the SOC effect, 14) which was later confirmed experimentally. 15) At higher pressures (∼18.9 Gpa), it was predicted that BiTeI changed from the P3m1 to P4/nmm phase and becomes superconducting by calculation, 16) after which the P4/nmm phase is experimentally proven to be a bulk superconductor, 17) showing the relevance of the structure to its unique properties. Further, topologically nontrivial centrosymmetric thin films of BiTeI 18) designed at 0 GPa exhibited various properties that were correlated to the structure.
In this work, we have explored the possible structures of BiTeI at 0 GPa and studied the electronic properties of several potential metastable structures of BiTeI via crystal structure prediction. These structures have high thermodynamic stability and show promise toward experimental synthesis. We have calculated band dispersion and a topological invariant (Z 2 ) of these structures and discovered two nontrivial topological structures. In addition to the pressure-induced topological phase 19) and centrosymmetric thin films, 20) BiTeI has other potential structures that can maintain a topological insulating state at 0 GPa and has a band gap of up to 860 meV. Through their spin texture, we confirmed the occurrence of Rashba splitting in the two structures and calculated their spin-orbit coupling parameters, and these structures may be conducive to the miniaturization of spintronic devices.
We computed 3000 crystal structures of BiTeI using global optimization and first-principles calculations at 0 GPa using the USPEX code. [21][22][23] The genetic algorithm of USPEX allows for an efficient global search for crystal structures without any prior knowledge. The simulation cells were 1-4 formula units of BiTeI. Because all structures had the same elemental ratio, the stability of the structures was evaluated by the enthalpy difference ΔH = H − H P3m1 between the predicted structure and the experimentally verified groundstate structure. 24) The energy evaluation and local geometric optimization of the structures were performed by the projector augmented plane-wave (PAW) methods 25) implemented in the Vienna Ab initio Simulation Package 26) (VASP). The Perdew-Burke-Ernzerhof (PBE) functional within the generalized gradient approximation (GGA), 27) and rev-vdW-DF2 28) which can offset the error caused by the dispersion force, were used to optimize the P3m1 structure that has been experimentally discovered. A plane-wave basis set is used and the kinetic-energy cutoff is set to 325 eV. The Brillouin zone of the P3m1 structure was sampled by 21 × 21 × 12 k-point meshes. The lattice parameters of the resulting optimized structure obtained by GGA (vdw-DF2) were a = 4.425(4.372) Å and c = 7.378(6.913) Å. Compared to the experimentally obtained lattice parameters of a = 4.339 Å and c = 6.854 Å, 24) the simulated ones were overestimated by 2% and 7.6% for the GGA-PBE and 0.7% and 0.8% for vdw-DF2. Therefore, the vdw-DF2 functional was used for all structural optimizations (see Sect. 1 of the supplementary material for more computational information on the other structures).
Here, the modified Becke-Johnson (mBJ) potential 29) was used to correct the band gap. The weight correlation constant was set to 1.1 and obtained an indirect band gap of 0.385 eV for the P3m1 structure, which is close to the experimental results of 0.38 eV 12) and 0.36 eV. 30) Thus, these parameters were used to correct the band gap and the band dispersion of all the structures.
We screened all the predicted structures to identify those with high stability and potential to be synthesized. Sun et al. 31) reported that 90% of all known inorganic metastable crystalline materials are 67 ± 2 meV/atom above the ground state. Herein, we used a more stringent criterion, i.e., structures with an enthalpy difference ΔH of less than 25 meV/atom from the ground state were considered to be highly stable, as shown in Table I.
Based on this criterion, five structures shown in Fig. 1 were identified finally (The cif files are given as the supplemental data). The experimentally verified P3m1 structure 12) [ Fig. 1(a)] was regarded as the ground state. In addition, we found the artificially designed two-dimensional topological insulator P3m1 20) [ Fig. 1

(b)] and confirmed its stability.
Three new structures were found, including two P1 structures [P1(I) and P1(II), ordered by their stability] and a layered structure P6 3 /mmc. Among these structures, P1(I) and P1(II) structures are non-centrosymmetric, whereas P6 3 /mmc is a centrosymmetric structure. These structures are characterized by the presence of corner-shared BiI 3 units that resemble regular triangular pyramids. The manner in which these units are arranged determines the overall structure, leading to the differences between those structures.
The P1(I) and P1(II) structures differ from the other structures as shown in Figs. 1(c) and (d). The distance between Bi atoms and other nearest-neighbor Bi atoms in all structures is about 4.3 Å, however, since there is only one atomic layer between the Bi layers of structures in the P1(I) and P1(II), the distances between different Bi layers are only 4.6 Å, while that of P3m1, P3 m1, and P6 3 /mmc are about 6.3 to 8.1 Å. P1(I) has a three-dimensional network, all the BiI 3 tetrahedrons are three-dimensionally connected. The P1(II) structure was obtained by swapping one of the Te atoms and I atoms within the P1(I) primary cell; the BiI 3 tetrahedrons are connected to each other in a plane.
The P6 3 /mmc structure [ Fig. 1(e)] showed inversion symmetry that originated from stacking in the c-axis. Similar to P3 m1, the two adjacent Te layers, and their corresponding BiTeI trilayers are stacked in the form of AB stacking. In AB stacking, wherein, atoms are displaced relative to the others by a vector equal to the hexagonal edge, giving rise to inversion symmetry in the structure. In contrast, the two adjacent I layers and their corresponding BiTeI trilayers in the P6 3 /mmc structure showed AA stacking, and the atoms have an identical horizontal coordinate. Figure 2 shows the band dispersions 33) of each structure by considering the SOC effect (For more high-symmetry k-paths defined by Setyawan et al. 34) which are often used, see Sect. 2 of the supplementary material). The band gap, largest spinorbit coupling parameter, and largest momentum offset are shown in Table II. The P3m1 structure shows spin splitting at the A-point as shown in Fig. 2(a), corroborated by the results of Ishizaka et al. 12) Since both the P1(I) and P1(II) structures are non-centrosymmetric; spin splitting was also present in their band dispersion as shown in Figs. 2(c) and 2(d). To shed light on the chemical bonding of the newly discovered structures, the effective charges of these structures were calculated by forces F under a finite electric field E through the Berry phase method 35) implemented in OpenMX. 36) The scalar charges, which are equivalent to one-third of the trace of the effective charge tensor ( ) were also calculated. One-third of the trace of the effective charge tensor is also called generalized atomic polar tensor charge, 37) which is useful for analyzing the polarity of populations. The results of the calculations To further verify the splitting type, PyProcar 38) was used to process and plot the spin texture of the P1(I) and P1(II) structures. The spin texture in the P1(I) structure around the Γ-point on the k x -k y plane is represented in Fig. 3(a), where x and y are in-plane coordinates for the a-b plane when the z-axis is set to the normal vector of the a-b plane, k x and k y are xand y-components of the wave vector, and the expectation values of the Pauli matrices σ x , σ y σ z at every k-point were evaluated. We investigated the spin texture at the energy level of E = E F + 0.65 eV above the degenerate bands at the Γ-point, where E F is the Fermi energy. The conduction band shows an oval internal and a nearly square external structure. The internal profile of the conduction band shows reversed spin polarization in relation to the external one, which is similar to the general characteristics of the Rashba spin splitting generated by spin Hamiltonian H so = a so (k x σ y − k y σ x ) where a so is spin-orbit coupling parameter. However, our analysis revealed that the spin texture of the P1(I) structure deviates from the typical Rashba splitting as it exhibits the out-ofplane components such as k x σ z and k y σ z and radial components such as k x σ x and k y σ y . The spin-orbit coupling parameter a so was estimated by using the expression a so = 2E so / k so , where E so and k so are the difference in energy and momentum, respectively, between the degenerate point of the time-reversal invariant momentum and its nearest-neighbor extreme point, as described in previous studies. 39,40) The calculations revealed that the spin-orbit coupling parameters of the Γ-point to Z and X directions were 0.42 eV·Å and 0.52 eV·Å, respectively.
For the P1(II) structure, splitting was found at both the Γand Z-points, and spin splitting around the Γ-point was shown in Fig. 3(b) In the same way as the P1(I) case, the spin as it has few out-of-plane components. This difference was explained through the analysis of the characteristics of crystal structures. The P1(II) structure has a structure close to that of R3m possessing the symmetry of C 3v : we confirmed that the P1(II) structure was classified into R3m with a tolerance of 0.07 Å through the FINDSYM. 41,42) From Fig. 2 in Sect. 4 of the supplementary material, we could see that there are three mirror planes and three-fold symmetry, C 3v . The previous work showed that the C 3v symmetry exhibits only the k x σ y − k y σ x term in the k-linear spin Hamiltonian. 43 Although the spin-orbit coupling parameters of the P1(II) structure were not as high as that of the P3m1 structure, the P1(II) structure has the largest offset k so of up to 0.077 Å −1 in the Z-U direction. According to the expression 44) for the differential phase shift Δθ = ΔkL, where Δk (=2k so ) is the differential wave vector between the up and down-spin electrons, and L is the distance traveled by electrons, as Δk increases, the distance required for the spin precession of the electron decreases. Therefore, a larger k so is usually considered to facilitate the miniaturization of the spin-type field effect transistor. Figures 2(b) and 2(e) show the band dispersion of the P3 m1 and P6 3 /mmc structures, respectively. The inversion symmetry in these structures leads to the disappearance of the spin splitting because the sum of potential gradients in each direction, which originated from the asymmetric distribution of internal charge, is zero. Both the P3 m1 and P6 3 /mmc structures showed a massive Dirac cone-like pattern at the Γ-point. The band gap of the P6 3 /mmc structure was 31meV, which was lower than that of the P3 m1 structure. Further, we studied the Z 2 topological invariant for the structures using DFT code OpenMX 45,46) and the post-process code Z2FH. 47,48) The Z 2 topological indices are summarized in Table II. By definition, the Z 2 topological invariant of the system can be determined by strong topological invariants n 0 , and n 1 , n 2 , n 3 , which are called weak topological invariants, using the following expression: Z 2 = (n 0 ; n 1 , n 2 , n 3 ). At 0 Gpa, both P3m1 structure and P6 3 /mmc structure belong to Z 2 = (1; 0, 0, 0), indicating that the P6 3 /mmc structure is a strong topological insulator, as same as the P3 m1 structure.
We also calculated the Z 2 invariants for the other structures and found that the P1(II) structure belongs to Z 2 = (0; 0, 0, 1), indicating that it is a weak topological insulator. Topological insulators are characterized by their unique high-mobility edge states that are topologically protected by time-reversal symmetry, which protects them from backscattering. This indicates that electrons can travel at high speeds along the surface of the material without dissipating energy due to scattering. However, since the band gap of known topological insulator materials is typically around 30-300 meV, 49) they can only function optimally at low temperatures. As the temperature increases, the band gap decreases and leads to internal conduction, severely limiting the practical applications of topological insulators. The band gap of the P1(II) structure being 860 meV is extremely rare, indicating that the device can function at high temperatures. This wide band gap can greatly expand the practical application scope of topological insulators.
Single crystal of BiTeI can be grown by the Bridgman method 16,50) and usually exists in the ground state; however, for the practical application of electronic devices, materials are often grown on a substrate. The structure of the crystal grown on a substrate is determined by the total free energy of the interface system. For general homophase and heterophase interface structures, the interface free energy, σ, can differ by order of magnitude, and the σ of heterophase and surface can be as high as 51) ∼124.8 meV·Å −2 . Therefore, using different substrate materials or substrates with thicknesses may result in the formation of metastable structures instead of groundstate structures. This phenomenon is known as "epitaxial stabilization" 52) and has been conversely applied to grow metastable materials such as oxide thin film 52) and superlattice material. 53) Except for epitaxial films, there are many other techniques to synthesize the desired structures, such as the recently reported synthesis of the metastable phase of bismuth selenide using quenching. 54) Compared to the results of previous studies, 18,20) we believe that growing the P3m1 structure on Bi 2 Te 3 substrates is a better choice than growing on BiTeI (P3m1) substrates. Lattice constants a of the P3m1 structure and P3m1 structure are 4.37278 Å and 4.38868 Å, respectively, while that of Bi 2 Te 3 is 4.386 Å, 55) which is closer to that of the P3m1 structure. Thus, the growth of the P3m1 and P3m1 structures on Bi 2 Te 3 induces a lattice strain to match the lattice constants, increasing a for the P3m1 structure by 0.3%, but decreasing that for P3 m1 only by 0.07%. Therefore, the free energy of the P3m1 structure grown on Bi 2 Te 3 substrates is 0.1 meV/atom lower than that of the P3m1 structure. In contrast, the free energy of the P63/mmc structure is 15.2 meV/atom higher than that of the P3m1 structure due owing to the AA stacking of the Bi-I layers. However, even for graphene, the free energy difference of AA stacking is as high as 12 meV/atom. 56) Therefore, the synthesis of the P6 3 /mmc structure may be practically realizable. For the P1 (I) and P1(II) structures, because their lattice parameters and structural properties differ, they could be grown on substrates with similar lattice parameters.
In conclusion, three novel BiTeI crystal structures were predicted by the genetic algorithm at 0 Gpa. The structural stability and electronic properties of the structures were investigated via first-principles calculations. Our systematic study showed that two of the structures were topological insulators at 0 GPa, including one strong topological insulator and one weak topological insulator with a wide band gap. Two of the structures had Rashba-type spin splitting, and one of them has a larger k so than the ground state (P3m1) structure. Our study can stimulate further investigation of potential metastable phases in other bismuth telluride halides BiTeX (X = Cl, Br) and promote the search for promising spintronic devices.