Importance of van der Waals interactions for ab initio studies of topological insulators
We investigate the lattice and electronic structure of the bulk and surface of the prototypical layered topological insulator BiSe using ab initio density functional methods. We show that inclusion of the van der Waals (vdW) interactions yields accurate values for both the lattice constants and distance between quintuple layers (basic structural units), in contrast to other standard approaches. For the topological surface state, inclusion of vdW interactions yields a Dirac cone with its tip just below the chemical potential, quasiparticle velocity very close to the experimental values, and a helical spin texture. Our results establish a framework for unified self-consistent first-principles calculations of topological insulators in bulk, slab and interface geometries, and provides the necessary first step towards ab initio modeling of topological heterostructures.
Most proposed applications of topological insulators (TIs) involve fabricating heterostructures combining these materials with other topologically trivial compounds, and creatively controlling the spin-momentum locked symmetry-protected states at the interface between the two. These interface states are often assumed to be helical and isotropically dispersing, in resemblance to their counterparts at surfaces. Hasan and Kane (2010); Qi and Zhang (2011) Theoretical models, however, show that symmetry breaking at the interface leads to substantial modifications of the properties of topological states. Zhang et al. (2012a); Asmar et al. (2017) In real systems such symmetry breaking interface potentials may originate from strain and stress due to lattice mismatch, broken bonds, buckling, and surface reconstruction; experiments indicate that they do indeed modify, sometimes drastically, the expected behavior of the topological states. Zhang et al. (2012b); Zeljkovic et al. (2015); Richardella et al. (2015) To develop a detailed understanding of the origin, form, and magnitude of these potentials, we first need a comprehensive picture of the structural and electronic properties of both the bulk and surface states which can be used as a reference and starting point for the interface problem.
Given the importance of this question, and noting that most topological materials are not strongly correlated, it would be reasonable to assume that the bulk, surface, and interface properties can all be understood from ab initio calculations. Surprisingly, however, few such studies accomplish a fully self-consistent picture. The reason for this is clear from reviewing existing treatments. Correct computation of the displacement of atoms in the vicinity of the interface requires, as the first step, accurate knowledge of the atomic positions and strain field in the bulk. For the prototypical topological insulator BiSe, an accurate self-consistent simultaneous determination of the bulk structure and electronic properties within density functional theory (DFT) approaches has been largely absent.
Here we fill this gap, and show that systematic inclusion of van der Waals (vdW) interaction into DFT-based approaches allows us to accurately and self-consistently obtain the lattice and electronic structure with correct and stable properties of the surface state. Our results demonstrate that accounting for vdW potentials in the ab-initio treatment of layered TIsFu and Kane (2007); Zhang et al. (2009) is a necessary first step towards a comprehensive understanding of such systems, which will lead to interface engineering and the design of topological electronic devices.
Ab initio approaches to BiSe and van der Waals forces.
BiSe, in addition to many other topological insulators, has a layered structure. The main structural unit is a “quintuple layer” (QL), which consists of three selenium and two bismuth layers arranged in the order SeBiSeBiSe, where the pairs of Se and Bi positions are related through inversion symmetry. The coupling between QLs is believed to be predominantly due to vdW forces. However, to our knowledge, this has not been properly implemented in DFT calculations for the bulk or surfaces. Surface calculations are typically done in the “slab” geometry, with a periodic arrangement of several QL-thick slabs separated by vacuum layers. Not using a sufficiently thick vacuum means that slabs interact with their periodic images through vacuum via long-range vdW interactions which affects the behavior of the surface states.
DFT studies make two choices at the outset. The first is the selection of the exchange correlation functional, with the Generalized Gradient Approximation (GGA) and the Local Density Approximation (LDA) most commonly used for BiSe. Methods such as the addition of the on-site Coulomb repulsion (LDA+U), use of hybrid functionals Crowley et al. (2015) that take into account a fraction of the exact Hartree-Fock exchange energy, and computing the electron self-energy with the renormalized Coulomb potential (GW) Yazyev et al. (2012); Nechaev et al. (2013), all aim to improve the accuracy of these approximations under various assumptions. The last of these is in principle the most precise, but both the GW method and its non-self-consistent counterpart, GW, are computationally expensive, rendering them impractical for surface and interface calculations Crowley et al. (2015). Therefore below we restrict ourselves to LDA, GGA, and LDA+U methods, which are sufficient for a weakly correlated material such as BiSe.
The second choice fixes the method for finding the optimal crystal structure. Many DFT studies of BiSe in bulk and at surfaces use experimentally determined lattice constants Yazyev et al. (2010, 2012); Crowley et al. (2015); Sun and Singh (2017), fixing the volume of the unit cell, and find the relaxed atomic positions within this cell under the constraint of respecting the crystal symmetry. The results reproduce salient features of the electronic spectra in bulk and at surfaces, although there is some debate about whether many body corrections are necessary to obtain the correct magnitude of the gap and character (direct vs indirect) Nechaev et al. (2013); Aguilera et al. (2013); Michiardi et al. (2014). In the same spirit, the study of interfaces has been largely limited to perfectly lattice matched cases Lee and Yazyev (2017), or to making a priori assumptions about structural changes at the boundary Zhang et al. (2016). Calculations are run with this method due to the fact that a structural optimization of bulk BiSe yields erroneous results: when a full geometry relaxation with LDA or GGA is run and the unit cell volume is allowed to change, the calculation yields lattice parameters different from those determined experimentally, implying that the strain field is not correctly determined using these methods Lind et al. (2005); Luo et al. (2012).
Previous attempts to include vdW interaction in the bulk Luo et al. (2012) and slab Zhang (2015); Govaerts et al. (2014) calculations have been inconsistent. Ref. Luo et al., 2012 performed structural optimization including vdW but omitting spin-orbit interaction (SOI), known to be critical for band formation and inversion in topological insulators; they then computed the electronic bands fixing this structure, employing GGA with SOI. Refs. Zhang, 2015; Govaerts et al., 2014 used slab geometries with small vacuum thickness, implying that the system was not free of electrostatic interaction of slabs with their images. As a consequence the obtained lattice constants and electronic dispersion differ from both experiment and our results presented below, and no bulk results are shown for comparison. We show that consistent inclusion of the vdW interaction in the GGA scheme improves agreement with experiment of the lattice constants computed with full geometric relaxation, by an order of magnitude, gives a reasonable value of the bulk band gap, and yields characteristics of the surface state very close to those observed experimentally.
We used a hexagonal unit cell for BiSe (see Fig. 1(a)) which is convenient for the description of the surface states. All calculations below were carried out using the Vienna Ab initio Simulation PackageKresse and Hafner (1993, 1994); Kresse and Furthmüller (1996a, b) (VASP), version 5.4.4. Crystallographic information is taken from experimental dataNakajima (1963) retrieved from Crystallography Open DatabaseMerkys et al. (2016); Gražulis et al. (2015, 2012, 2009); Downs and Hall-Wallace (2003). We used Project Augmented Wave (PAW) potentialsBlochl (1994); Kresse and Joubert (1999) for Bi () and Se (), for a total of 48 electrons, and a plane-wave basis. Convergence tests revealed that a -centered k-point grid of 11x11x11 k-points and an energy cut-off of 450 eV for the plane wave basis are sufficient for high accuracy results. We performed fully relativistic calculations which include SOI. The convergence threshold for energy is taken to be eV. Band structures are plotted with data processed using vaspkitWang (2013). We used GGA–PBEPerdew et al. (1996, 1997), LDAPerdew and Zunger (1981), LDA+U and GGA+vdW, with full relaxation. We compared the results of two methods including van der Waals interactions: semi-empirical DFT–D2 Grimme (2006) and DFT–D3 Grimme et al. (2010) methods. Semi-empirical methods add to the total energy corrections proportional to (DFT–D2) with additional terms varying as (DFT–D3) for each pair of atoms which are separated by less than the cutoff distance. We used the default cutoff distance of 50 Å, and checked that cutoff radii of 40Å and 65Å do not change the results. To reduce the contribution from pairs of atoms that are bonded covalently these methods use a short-distance damping function, and we used the original Fermi type function Grimme (2006) for DFT–D2, and the D3(zero) form for the DFT-D3 calculation.
Our results for the structural optimization are shown in Table 2. As expected, LDA overbinds the electrons leading to a 1.3% contraction of the lattice constants compared to their experimental value. This reduction is in quantitative agreement with the 4% volume change of the unit cell found in Ref. Lind et al., 2005, albeit the values for the lattice constants differ due to different optimization procedures. The unit cell volume is overestimated with GGA by about 7% in our calculation vs. almost 10% in Ref. Lind et al., 2005. Ref. Luo et al., 2012 found the cell volume overestimated even more using PBE, but that is mostly due to a significant elongation of the -axis lattice constant, perhaps related to the effectively two-dimensional k-point mesh (13x13x1) used in that work. Both GGA and LDA give large errors for the bulk energy gap, as is common for small gap semiconductors with strong spin-orbit coupling.
Notably, using both GGA and LDA we find that the difference between the experimental values and those obtained from first principles is much greater for the distance between adjacent QLs (see Fig. 1(a)), in Table 2, than for the lattice constants. Since the QLs represent a system with an effectively closed shell Lind et al. (2005), the bonding between the QLs is due to van der Waals forces, and this result immediately suggests the critical role of the vdW interaction for structural optimization.
To become convinced that intra-QL Coulomb correlations cannot account for the discrepancy, we performed LDA+U calculations with the on-site repulsion on Bi orbitals, see Table 2. For , the lattice constants and are very close to their experimental values, but the inter-QL distance still significantly deviates from that in experiment, and is only improved for , where the lattice constants increase. In contrast, the gap magnitude is closest to experiment for . Therefore no single value of consistently improves the results. Moreover, the observed trends imply that the agreement with experiment is accidental: increasing values of tend to collapse the energy gap, in contrast to the physical expectation that Coulomb repulsion localizes corresponding orbitals and pushes bands apart. The situation is even worse for GGA+U as we very quickly reach gap collapse and metallicity. We therefore conclude that the error in the interlayer distance must be related to the long-range part of the interaction.
Including van der Waals in the structural optimization using the Grimme scheme Grimme (2006), we find a dramatic improvement (an order of magnitude) in the agreement with experiment for values of the lattice constants, and the interlayer spacing, see last column of Table 2. The inter-QL spacing deviates from experiment by only , while the lattice constants match the experimental values almost exactly. The obtained band gap (indirect) is much closer to the experimental value than that obtained using other functionals, see Table 2. The band structure shown in Fig. 1(c) exhibits, albeit not very strongly, a characteristic ‘camelback’ feature in the valence band at the center of the Brillouin Zone () due to Spin-Orbit Coupling. This feature is sensitive to the choice of the approximation: it nearly vanishes in GGA, and is over-emphasized in LDA to the extent that the conduction band acquires this feature as well. It has been argued that the many body GW corrections remove this feature and produce a direct gap that agrees with experiment Michiardi et al. (2014), but resolving this controversy is not the main focus of our work. We note, however, that the inclusion of a GW correction into LDA that “straightens” the camelback feature also leads to a substantial reduction of the magnitude of the energy gap Yazyev et al. (2010); Aguilera et al. (2013); Michiardi et al. (2014) (overestimated in LDA by over 50%, see Table 2), while in our case the gap is much closer to the experimental value. Note also that the band structure exhibits substantial energy dispersion along the direction, which suggests that the k-point mesh with only one point along used in some previous DFT-vdW calculations Luo et al. (2012); Zhang (2015) is insufficient for accurate description of BiSe.
In summary, we showed that the addition of vdW corrections compensates for the severe error in the interlayer distance in the standard approaches, and yields correct structural parameters as well as improved band gap and electronic structure without the computational cost of many-body approaches. We conclude therefore that inclusion of vdW forces is essential in any full relaxation calculations of the layered TI structures.
Slab calculations and the surface states.
To determine the structure of the surface states we used the same GGA+vdW method with the D2 choice for the vdW correction. We found that for 5QL-thick ( 50 Å) and thicker slabs the surface states at the opposite faces hybridize sufficiently weakly so that there is no observable gap in the Dirac spectrum. We also determined that to model surface states we need to include an amount of vacuum that is roughly equal to twice the slab thickness in order to avoid electrostatic interaction of the slab with its own periodic images, which notably also generates a gap in the spectrum of the surface states. In our calculations the outermost QLs of the slab (five atomic layers closest to each surface) are allowed to relax, while the atoms in the “bulk” part of the slab are kept fixed. We find that the atomic displacements are of order mÅ, negligible for the electronic structure calculations presented below.
The slab band structure is shown in Figure 2(a), and shows quasi-linearly dispersing Dirac-like states in the bulk gap arising from the surface states. The vertex of the Dirac cone lies just below the chemical potential. This is in qualitative agreement with experiment Zhang et al. (2009), although self-doping effects likely affect the location of the Dirac point. Importantly, we find the velocity of the Dirac quasiparticles to be m/s, which is much closer to the experimentally determined values of m/s Devidas et al. (2017); Pertsova et al. (2016); Kuroda et al. (2010); Chang et al. (2015); Cao et al. (2013) than the values obtained in other ab initio calculationsYazyev et al. (2010); Crowley et al. (2015).
We find that the surface states are nearly perfectly helical, as shown in Fig. 2(b). The spins are predominantly normal to the direction of the momenta, and reverse direction upon crossing the Dirac point. The insets show the calculated spin orientation along two different constant energy surfaces. We find that the spins lie in the plane with negligible out-of-plane component. We also find no evidence of hexagonal warping of the constant energy contours at higher energies (for example, at eV), in agreement with experiments on BiSeNomura et al. (2014); Kuroda et al. (2010).
The spatial profile of the surface states is shown in Fig. 3. The decay length of the surface state is of order of a single QL, and, for our 7QL-thick slab, the weight essentially vanishes at the center. For comparison, we show the amplitude of the wave function for one of the bulk band states, which is reduced near the surface, and extends throughout the bulk.
Summary. We performed a comprehensive high precision ab-initio analysis of the lattice structure and electronic properties of the prototypical topological insulator BiSe comparing different approaches such as LDA and GGA-PBE density functional methods with and without van der Waals interactions. Our results show that inclusion of van der Waals interactions is critical for ab initio studies of layered topological materials. The key observation is that the distance between the closed shell-like layers determined from standard DFT methods without vdW corrections yields a large deviation from the experimental values. Using LDA+U methods produces inconsistent results with non-physical trends. Inclusion of the vdW interaction improves the agreement with experimental structural constants for the bulk material by nearly an order of magnitude, significantly improves the value for the energy gap, and yields dispersion and spin structure of the topological surface state that closely matches experimental results. Our results establish a pathway for reliable self-consistent determination of the surface and interface properties of topological materials, and for ab initio analysis of prototype topological devices including stress, strain, and symmetry breaking effects.
Acknowledgments. This research was supported by NSF via Grant No. DMR-1410741 (K. S. and I. V.) and by the U.S. Department of Energy under EPSCoR Grant No. DE-SC0012432 with additional support from the Louisiana Board of Regents (W. A. S.).
- Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Reviews of Modern Physics 82, 3045 (2010).
- Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Reviews of Modern Physics 83, 1057 (2011).
- Zhang et al. (2012a) F. Zhang, C. L. Kane, and E. J. Mele, Phys. Rev. B 86, 081303(R) (2012a).
- Asmar et al. (2017) M. M. Asmar, D. E. Sheehy, and I. Vekhter, Phys. Rev. B 95, 241115(R) (2017).
- Zhang et al. (2012b) Q. Zhang, Z. Zhang, Z. Zhu, U. Schwingenschlögl, and Y. Cui, ACS Nano 6, 2345 (2012b), .
- Zeljkovic et al. (2015) I. Zeljkovic, D. Walkup, B. A. Assaf, K. L. Scipioni, R. Sankar, F. Chou, and V. Madhavan, Nat. Nanotechnol. 10, 849 (2015).
- Richardella et al. (2015) A. Richardella, A. Kandala, J. S. Lee, and N. Samarth, APL Mater. 3, 083303 (2015).
- Fu and Kane (2007) L. Fu and C. L. Kane, Physical Review B 76, 045302 (2007).
- Zhang et al. (2009) H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Nature Physics 5, 438 (2009).
- Crowley et al. (2015) J. M. Crowley, J. Tahir-Kheli, and W. A. Goddard, The Journal of Physical Chemistry Letters 6, 3792 (2015).
- Yazyev et al. (2012) O. V. Yazyev, E. Kioupakis, J. E. Moore, and S. G. Louie, Physical Review B 85, 161101(R) (2012).
- Nechaev et al. (2013) I. A. Nechaev, R. C. Hatch, M. Bianchi, D. Guan, C. Friedrich, I. Aguilera, J. L. Mi, B. B. Iversen, S. Blügel, P. Hofmann, and E. V. Chulkov, Physical Review B 87, 121111(R) (2013).
- Yazyev et al. (2010) O. V. Yazyev, J. E. Moore, and S. G. Louie, Physical Review Letters 105, 266806 (2010) .
- Sun and Singh (2017) J. Sun and D. J. Singh, Journal of Applied Physics 121, 064301 (2017) .
- Aguilera et al. (2013) I. Aguilera, C. Friedrich, G. Bihlmayer, and S. Blügel, Physical Review B 88, 045206 (2013).
- Michiardi et al. (2014) M. Michiardi, I. Aguilera, M. Bianchi, V. E. de Carvalho, L. O. Ladeira, N. G. Teixeira, E. A. Soares, C. Friedrich, S. Blügel, and P. Hofmann, Physical Review B 90, 075105 (2014).
- Lee and Yazyev (2017) H. Lee and O. V. Yazyev, Phys. Rev. B 95, 085304 (2017).
- Zhang et al. (2016) J. Zhang, J. P. Velev, X. Dang, and E. Y. Tsymbal, Physical Review B 94, 014435 (2016).
- Lind et al. (2005) H. Lind, S. Lidin, and U. Häussermann, Physical Review B 72, 184101 (2005).
- Luo et al. (2012) X. Luo, M. B. Sullivan, and S. Y. Quek, Physical Review B 86, 184111 (2012).
- Zhang (2015) Y. N. Zhang, The Journal of Chemical Physics 143, 151101 (2015) .
- Govaerts et al. (2014) K. Govaerts, K. Park, C. De Beule, B. Partoens, and D. Lamoen, Physical Review B 90, 155124 (2014).
- Nakajima (1963) S. Nakajima, Journal of Physics and Chemistry of Solids 24, 479 (1963).
- Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558(R) (1993).
- Kresse and Hafner (1994) G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994).
- Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Comput. Mat. Sci. 6, 15 (1996a).
- Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
- Merkys et al. (2016) A. Merkys, A. Vaitkus, J. Butkus, M. Okulič-Kazarinas, V. Kairys, and S. Gražulis, J. Appl. Cryst. 49 (2016).
- Gražulis et al. (2015) S. Gražulis, A. Merkys, A. Vaitkus, and M. Okulič-Kazarinas, J. Appl. Cryst. 48, 85 (2015).
- Gražulis et al. (2012) S. Gražulis, A. Daškevič, A. Merkys, D. Chateigner, L. Lutterotti, M. Quirós, N. R. Serebryanaya, P. Moeck, R. T. Downs, and A. Le Bail, Nucleic Acids Research 40, D420 (2012).
- Gražulis et al. (2009) S. Gražulis, D. Chateigner, R. T. Downs, A. T. Yokochi, M. Quirós, L. Lutterotti, E. Manakova, J. Butkus, P. Moeck, and A. Le Bail, J. Appl. Cryst. 42, 726 (2009).
- Downs and Hall-Wallace (2003) R. T. Downs and M. Hall-Wallace, American Mineralogist 88, 247 (2003).
- Blochl (1994) P. E. Blochl, Phys. Rev. B 50, 17953 (1994).
- Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
- Wang (2013) V. Wang, “Vaspkit — postprocessing tool for vasp code,” (2013).
- Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- Perdew et al. (1997) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 78, 1396 (1997).
- Perdew and Zunger (1981) J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
- Grimme (2006) S. Grimme, J. Comp. Chem. 27, 1787 (2006).
- Grimme et al. (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, The Journal of Chemical Physics 132, 154104 (2010).
- Devidas et al. (2017) T. R. Devidas, E. P. Amaladass, S. Sharma, A. Mani, R. Rajaraman, C. S. Sundar, and A. Bharathi, Materials Research Express 4, 026101 (2017).
- Pertsova et al. (2016) A. Pertsova, C. M. Canali, and A. H. MacDonald, Physical Review B 94, 121409(R) (2016).
- Kuroda et al. (2010) K. Kuroda, M. Arita, K. Miyamoto, M. Ye, J. Jiang, A. Kimura, E. E. Krasovskii, E. V. Chulkov, H. Iwasawa, T. Okuda, K. Shimada, Y. Ueda, H. Namatame, and M. Taniguchi, Physical Review Letters 105, 076802 (2010).
- Chang et al. (2015) C.-Z. Chang, P. Tang, X. Feng, K. Li, X.-C. Ma, W. Duan, K. He, and Q.-K. Xue, Physical Review Letters 115, 136801 (2015).
- Cao et al. (2013) H. Cao, S. Xu, I. Miotkowski, J. Tian, D. Pandey, M. Z. Hasan, and Y. P. Chen, physica status solidi (RRL) - Rapid Research Letters 7, 133 (2013).
- Nomura et al. (2014) M. Nomura, S. Souma, A. Takayama, T. Sato, T. Takahashi, K. Eto, K. Segawa, and Y. Ando, Physical Review B 89, 045134 (2014).