Ab-initio transport across Bismuth Selenide surface barriers
We investigate the effect of potential barriers in the form of step edges on the scattering properties of BiSe(111) topological surface states by means of large-scale ab-initio transport simulations. Our results demonstrate the suppression of perfect backscattering, while all other scattering processes, which do not entail a complete spin and momentum reversal, are allowed. Furthermore, we find that the spin of the surface state develops an out of plane component as it traverses the barrier. Our calculations reveal the existence of quasi-bound states in the vicinity of the surface barriers, which appear in the form of an enhanced density of states in the energy window corresponding to the topological state. For double barriers we demonstrate the formation of quantum well states. To complement our first-principles results we construct a two-dimensional low-energy effective model and illustrate its shortcomings. Our findings are discussed in the context of a number of recent experimental works.
Bismuth selenide has emerged as the prototypical topological insulator material due to a single Dirac cone in the surface band structure and a relatively large bulk band gap. In 2009, concurrent theoretical (1) and experimental (2) works revealed the topological insulator phase of BiSe. Since then many fundamental properties of topological states have been demonstrated in this material, which has been called the Hydrogen atom of topological insulators. (3); (4)
In recent years, there has been a rapid expansion in the number of scanning tunneling microscope (STM) experiments on the BiSe and the closely related BiTe surface. Impurities on bismuth selenide have been imaged and scattering mediated by bulk states has been observed. (5); (6); (7); (8) Additionally there have been studies of dopants on the bismuth telluride surface (9); (10) and, interestingly, a bound state at a surface step of BiTe has also been found. (11) On the theoretical front, there have been several efforts towards modeling scattering of these surface states from perturbation theory by employing Dirac-like model Hamiltonians and by imposing symmetry considerations. (12); (13); (14) Furthermore, a study of robustness of surface states against on-site disorder by employing first-principles calculations was also reported. (15) The problem of scattering at a monolayer-bilayer graphene junction has also been investigated. (16)
In this paper, we investigate the effect of step barriers at the BiSe(111) surface on the scattering properties of the topological states by means of ab-initio transport calculations. We find that, due to the spin-polarized helical nature of the surface band, there is no scattering for normal incidence, since a reflection would entail a backscattering. However, as one moves to non-normal incidence, scattering is revealed. This is because the spins of the counter-propagating channels are no longer anti-parallel. An analysis of the local density of states reveals that the surface barrier strongly affects the spin of the surface state, in particular allowing an out of plane spin component, which is negligible in the absence of the barrier. In order to compare to our ab-initio results we have constructed a potential barrier model based on a simple Dirac Hamiltonian for the surface states. This is solved for barriers of various shapes and a comparison is made with our first-principles calculations. We note in passing that, although our ab-initio calculations have been performed for the particular case of bismuth selenide, we expect the same qualitative results to also hold for step edges perpendicular to directions without hexagonal warping in BiTe and for other related materials like BiTeS and TlBiSe.
Following this introduction, the remaining of the paper is organized as follows. We begin by describing our computational methods, in particular we outline the procedure for performing the transport calculations. We then study scattering originating from a single surface barrier by analyzing the transmission and the densities of states. Intriguingly our calculations reveal a bound state in the vicinity of the barrier, and we study its energy dispersion. From there, we move on to construct a low-energy model for the scattering problem, with barrier strength and width extracted from our ab-initio results. Next we look at the analogous problem in the presence of double barriers of different lengths. Notably, we find an energy splitting of the bound state when the states at the two barriers interact directly as in the case of a short double barrier, as well as when the bound states couple with quantum well states formed in the case of a longer double barrier. Finally we summarize our results and conclude.
Ii Computational Methods
Our transport calculations have been performed by using the Smeagol code, which combines the density functional theory (DFT) numerical implementation contained in the siesta code (17) with the non-equilibrium Green’s function method for electron transport. Here we briefly outline the calculation procedure and refer the readers to References [(18); (19); (20)] for a more detailed exposition. In Smeagol semi-infinite electrodes are attached to a central scattering region by means of self-energies. The calculation of the self-energies of the leads is performed by using a singular value decomposition-based, robust and efficient algorithm, which overcomes the problems related to recursive methods. (20) The Hamiltonian needed for the algorithm is calculated by using an equivalent infinite bulk system. The transport calculation proceeds by using the density matrix and the Kohn-Sham Hamiltonian obtained from siesta. The self-energies for the leads are added to the Hamiltonian of the scattering region, , and the non-equilibrium Green’s function is obtained by direct inversion
where is the self-energy of the left-hand side (=L) and right-hand side (=R) lead. The charge density is then calculated by integrating the non-equilibrium Green’s function along a contour in the complex energy plane
where is the lesser Green’s function for the transport problem. The Fermi functions for the leads are denoted by and are the broadening matrices. In order to perform the contour integral in Eq. (2) we use 16 energy points in the complex semicircle, 16 points along the line parallel to the real axis and 16 poles. The density matrix calculated in Eq. (2) is used by siesta to re-evaluate the Kohn-Sham Hamiltonian, and such a procedure is iterated until self-consistency is obtained. Once convergence is achieved, the relevant quantities like the transmission function, , and the density of state (DOS), , are calculated,
where stands for the trace, is the spectral function and is the overlap matrix.
In all calculations spin-orbit interaction is included by means of an on-site approximation (21) and the Perdew-Burke-Ernzerhof generalized gradient approximation (GGA) to the exchange-correlation functional is employed. (22) We have used a double- polarized basis set and a real space mesh cutoff of 300 Ryd. For slab calculations a minimum of 25 Å vacuum region has been included to prevent spurious interaction between periodic replicas. We use -point mesh to obtain the self-consistent potential (here is the direction perpendicular to the transport direction in the plane of the slab, is along the slab height and is the transport direction). When calculating the integrated transmission and DOS we use -points along the direction. Periodic boundary conditions have been considered in the plane orthogonal to the transport direction, while using open boundary conditions along the transport direction allows us to simulate a single scatterer, which, in this particular case is a surface step.
The unit cell used for the leads is shown in Fig. 1(a). We note that by adding the self-energies obtained from such quasi two dimensional leads to the Hamiltonian of the scattering region one avoids any finite size effects. It consists of a three quintuple-layers (3-QL) thick slab terminated on both sides by Se atoms, as found experimentally. For the slab we use the experimental lattice constants. The corresponding band structure is shown in Fig. 1(b). Note that there is band folding as a consequence of the doubling of the BiSe primitive unit cell. The bandstructure reveals the Dirac cone and the helical states consistent with earlier studies. (1) It should also be noted that there is a small but finite gap (of the order of 0.015 eV) at the point in the cone due to interaction between the two surfaces at opposite sides of the slab. However, this small gap does not affect our analysis of the topological states at higher energies, since the tunneling between the two surfaces is negligible as we elaborate in the next section. The transport setup for single and double barrier is shown in Fig. 1(c) and (d) respectively. We consider a single QL-high barrier on a 3-QL thick slab. The step edge is extended along the direction and the transport is along the orthogonal direction of the primitive BiSe Brillouin zone. The scattering region has a length of 198.87 Å. For the single barrier case the 4-QL region extends over about half the length of the scattering region. For the double barrier setup we investigate two barrier lengths, where the step extends over a region of 49.72 Å in the shorter case and is 149.16 Å for the longer one. Note that such large cells comprising of a few thousand atoms require an accurate order- algorithm as the one available in Smeagol. (23)
Iii Scattering from a single barrier
We begin our analysis by looking at the transport across a single surface barrier [see Fig. 1(c)], for which the transmission function is shown in Fig. 2(a) as a function of energy and for different values of the component of the wave-vector. At normal incidence (), the surface states are perfectly transmitted, , due to their helicity. As such, our first-principles calculations confirm Klein tunneling. (24) The transmission of bulk states, however, is reduced by the presence of the step edge. In contrast, as soon as one moves away from normal incidence, the transmission is no longer integer-valued. In particular it dips below , indicating substantial scattering. Note that the drop in transmission at at eV is merely due to the small gap in the band structure due to the finite thickness of the slab. Fig. 2(b) shows the total transmission obtained by integrating over all angles of incidence, namely , where is the length of the Brillouin zone. Notably retains the characteristic “V-shape” associated with the linear Dirac cone-like bands, despite the presence of the barrier. Overall we can conclude that the total transmission in presence of the barrier is quite close to the one for the unperturbed slab [compare the red and black curves in Fig. 2(c)]. For comparison, we have also performed calculations for steps running along the direction (with transport along ). Since the hexagonal warping effect, particularly at energies close to the Dirac crossing, is quite small in BiSe, we find results, which are very similar to the ones obtained for steps along the direction. Hence, in the rest of this paper we focus our attention on the latter.
At non-normal incidence the spin projections of the surface states counter-propagating at a given edge are no longer anti-parallel and thus backscattering becomes allowed, even in the absence of a perturbation that breaks time-reversal symmetry. We note that, although spin-orbit coupling mixes the spin components, one can still define spin components along different directions by using a projection onto the three Pauli matrices and the identity matrix . The situation is schematically illustrated in Fig. 2(d), and its consequences are demonstrated in Fig. 2(c), where we plot the transmission across the surface barrier as a function of at different energies. Clearly is reduced as increases, which is expected from argument related to the spin projections of the two counter-propagating surface states. At larger incidence angles the transmission tends towards the residual value of one, since a perfectly transmitted surface state is present at the opposite side of the slab (no scattering center is present on the opposite surface). If one increases even further, the band edge for the Dirac cone at both surfaces is reached, and the transmission abruptly goes to zero. It can be shown that the maximum scattering amplitude is proportional to , where is the angle between the spin directions of the counter-propagating surface states. (25) Note that at higher energies, the transmission persists at values around the unperturbed one, , for larger incidence angles. This is because as one moves the Fermi level at higher energy, the Fermi circle gets larger. Consequently, the same corresponds to a smaller incidence angle.
In STM experiments, one measures the oscillations in the electron density in order to study the scattering arising from surface modifications, for example from surface steps as studied in Ref. (26). A Fourier transform of the density yields the characteristic frequencies of its oscillations, i.e, gives the scattering wavevectors, ( and are incident and reflected wavevectors, respectively). In Fig. 3 we plot the density of states projected (PDOS) onto the surface atoms along the scattering region. At no oscillations in PDOS are seen after reflection from the step edge. However, moving away from normal incidence, the above-mentioned oscillations begin to appear. We remind the reader that along the transport direction () the use of self-energies corresponding to the left-hand side and right-hand side electrodes makes the system infinite but non-periodic. Thus no finite size effects, orthogonal to the step direction, are seen on the density oscillations. The scattering vectors can be obtained by performing a Fourier transform of the DOS along the long flat region adjacent to the barrier. At , expectedly there are no prominent scattering processes. As one moves to Å, there appears a dominant scattering wave-vector in the Fourier transform starting at 0.1 eV and extending upwards in energy, as shown in Fig. 3(b). This corresponds to backscattering at a non-normal incidence angle. Furthermore, this can be mapped to band structure along the transport direction, where a band starting at the same energy is present. The average over , however, reveals no scattering on this scale, even though there is a clear back scattering at individual . In order to accurately resolve the small density oscillations above the average, one would need to consider many more -points in the calculation. This is computationally prohibitively expensive for the system sizes considered here, and a more detailed investigation of the oscillations will be reported elsewhere. (27) For all three cases we also plot the transmission as a function of energy, for comparison with the surface PDOS.
Fig. 3 also makes apparent the band bending (of the order of 0.04 eV) introduced by the step. We will show in the next section that such band bending close to the step is a crucial ingredient for constructing a scattering model. Far enough from the step, however, the PDOS reverts to the unperturbed value within 40 Å, consistent with experimental observation. (28)
In contrast to similar steps on the Sb(111) surface, (29); (30) in BiSe we find bound states close to the step edge and penetrating into the barrier (with an exponentially damped oscillating amplitude). These exist over the entire energy window in which the surface states are present. Similar features with enhanced DOS have been measured by Alpichshev et al. (11) around surface barrier at the BiTe surface. Importantly such a bound state was not ascribed to the warped band structure of BiTe. Our results point towards a similar bound state in BiSe as well. In the experiment, no information could be obtained about the DOS on the lower side of the step. Our calculations in fact reveal that the state exists only on the higher side of the barrier, and the lower side has no such features. We have also calculated the energy dispersion of this state along the direction perpendicular to the transport. We plot the energy and dependence of the PDOS on the Se atom at the barrier [shown in Fig. 4(b)] and on a surface atom 50 Å away from the barrier [Fig. 4(c)], and compare them to the band structure for the perfect periodic system [Fig. 4(a)]. For the atom present at the barrier we find additional pair of states outside the unperturbed Dirac bands, which is consistent with topological band theory. These additional states merge into the Dirac cone at eV and produce an enhanced PDOS around that energy. Away from the barrier, however, the PDOS is very similar to that of the unperturbed system, consistent with the pair of bound states being present only close to the barrier. We believe that these predictions of the bound state in BiSe and its energy dispersion may find verification in future experiments.
By analyzing the PDOS of the atoms at the bottom surface we have checked that significant scattering occurs only at the top one, i.e., it is caused by the presence of the step edge and not due to the tunneling back to the bottom surface. In Fig. 5 we plot the PDOS on the atoms present at the bottom surface at normal incidence and at a representative value of Å for oblique incidence. For both cases we can see absence of density oscillations in the energy range corresponding to the surface bands. Notably no signature of the bound state is also observed. Furthermore, we evaluate the Fourier transform of the PDOS in the flat region next to the step and find no features which may be mapped back to the scattering wavevector , in contrast to the case of the top surface.
The local density of states (LDOS) associated to electronic states incoming from the left-hand side lead at 0.175 eV above the Fermi level are shown in Fig. 6. (31) These clearly illustrate the three-dimensional nature of the path that electrons must traverse while crossing the barrier. The spin projections of the LDOS at and Å are shown in the left and right panels, respectively. In contrast to pristine bismuth selenide the spins of the helical surface states are no longer confined to the plane of the slab. In the vicinity of the barrier they rotate out of the plane (the component becomes finite). The LDOS at the bottom unperturbed surface provide a convenient comparison to the pristine surface, albeit with the spin directions reversed. At Å, the and components are dominant for the bottom surface, while the step edge introduces a component along the direction comparable with the other two, at the top surface. A zoom close to the step shows a large DOS close to the step edge, which is due to the bound state.
Iv A low-energy model
In order to interpret our ab-initio results we construct a simple potential barrier model for the scattering problem. The surface states are described by a Dirac Hamiltonian (1)
where the potential profile is shown in Fig. 7(a). The values of eV and eVÅ, are obtained from our first-principles band structure. Here we consider only the upper part of the cone, i.e., . The corresponding eigenstate is given by,
One can then use the wave-function continuity conditions at the potential steps to solve for the transmission and reflection coefficients in a straightforward manner. The potentials in the 4-QL and 3-QL leads, respectively and , are nearly identical and are set to zero. is the potential associated to the barrier and extends over a length , while is the band bending, which is finite over a distance . The calculated transmission curves are plotted in Fig. 7(b) for eV and in Fig. 7(c) for . The shape of the transmission function is much closer to that obtained from the ab-initio calculations for finite (this value of is chosen from our first principles results), as compared to the situation where . While this comparison does not provide definite evidence of importance of band bending, it serves as an illustration that it is one of the factors which need to be considered while performing a quantitative modeling of step edges on topological insulator surfaces. Although it appears that this simplified model can qualitatively reproduce the transmission obtained from first-principles, a more careful analysis shows that it neglects a number of important aspects of the scattering problem. It does not take into account the three-dimensional nature of the barrier, and as a consequence it cannot capture the change in spin orientation of the surface states near the barrier. Moreover it needs as an input the values of the scattering potentials, which an atomistic description is capable of providing, while also capturing the fine details of the scattering process. We also note that several models have been proposed to study topological states on a curved surface. These predict no backscattering at any angle from hyperbolic steps. (32); (33) Unfortunately these models are not valid for atomic-scale abrupt steps that we have studied in this work.
V Scattering from double barriers
We now analyze the scattering properties of double barrier structures constructed over the BiSe(111) surface. The scattering region is shown in Fig. 1(d) for the shorter surface barrier. This time the scattering structure is connected on both sides to two identical semi-infinite leads (3-QL slabs). As before, we begin by looking at the transmission across the surface as shown in Fig. 8(a). Again counter-propagating spin-momentum-locked states yield a perfect transmission at normal incidence. As discussed for the single barrier case, at finite the transmission is then reduced. However, in contrast to the previous analysis, there are also resonant energies at which the transmission reaches up the value of two, i.e., there is no reflection. At these particular energies the system displays Fabry-Perot resonances, which are characteristic of one-dimensional scattering from double potential barriers. In Fig. 8(b) we plot the transmission as a function of the incident for different energies. Away from the resonances the transmission shows again a cosine-like behavior with transmission going down to as the incidence angle increases ( gets larger). At even larger (not shown) the transmission drops down to zero when the band edge for the Dirac cone is reached at the bottom surface, similar to the case of single barrier.
The -resolved and total DOS projected on the surface atoms is plotted in Fig. 9, where the bound state can be clearly seen in the 4-QL region extending from 10 Å to 60 Å. The DOS associated to such bound state oscillates and decays towards the center of the quantum well defined by the two barriers at the step edges. A band bending similar to that observed for the single barrier is also seen for this particular barrier configuration. The interaction between the bound states localised at the two barriers splits them in energy, creating alternating high and low DOS as one moves up along the energy axis. Another noticeable feature is a state localized in the 4-QL region at around 0.1 eV [see Fig. 9(b)]. This is an additional state in the 4-QL slab, which is decoupled from the 3-QL leads. The same state is absent in the case of a single barrier produced by a step edge between a 3-QL and a 4-QL semi-infinite lead. The Fourier transforms of the DOS display similar features as those shown in Fig. 3. However, in the double barrier case the resolution is improved over that of the single barrier structure since we now have more atoms along the flat region next to the barrier.
We also study the energy dispersion of the quasi-bound state obtained at the barrier, by calculating the PDOS on the Se atom at the barrier, as a function of and momentum along the step (). This is shown in Fig. 4(d), with a comparison to the band structure of the unperturbed periodic system. Apart from the Dirac bands, additional states, dispersing along are visible at the interface. These have a dispersion very similar to the case of a single barrier [see Fig. 4(b)]. However, some additional features are seen when this pair of states mixes with the Dirac bands, with an alternating pattern of higher and lower PDOS being visible. This is due to the interaction between the bound states at the two barrier edges. Away from the interface the PDOS and the dispersion reverts to that of the pristine system with only the Dirac bands being present.
Note that for this particular chosen length of the double barrier there are no quantum well states formed inside the 4-QL region. However, for a longer barrier the quantum well states appear, as demonstrated by the PDOS on the surface atoms at two different for a barrier of length 149.16 Å (see Fig. 10). At normal incidence no quantum well states can be formed in the energy window of the surface state, since the two surface states have opposite spin projections leading to no interference. In contrast, at finite quantum well states appear (e.g. a nodeless state at around 0.13 eV and a single-node state at around 0.16 eV). However, the behavior of these states near the edges of the barrier is different from usual because of the presence of the bound state. In fact, these quantum well states interact with the bound states at the edges of the barrier resulting in an energy splitting of the bound state. We observe splitting of the bound states in both the short and the long double barrier, in the former case due to the interaction between the bound states located at the two edges of the 4-QL region, while in the latter due to the bound state interacting with the quantum well state within the barrier.
We have used ab-initio transport theory to study scattering to both single and double barriers of the topological protected states present on a BiSe(111) surface. In particular we have studied the dependence of the transmission on the angle of incidence and the electron energy. At normal incidence our first principles approach confirms Klein tunneling. Furthermore, we have calculated the density of states projected on the surface atoms and found bound states localised only on the higher side of the barrier. Thus our local density of states plots make apparent the three-dimensional nature of the scattering problem, in which the spins of the surface states are no longer confined to the plane of the topological insulator slab. We have also constructed a simplified potential barrier model using linear Dirac bands to compare with our first principles calculations. Throughout the paper we have placed our results in the context of recent experimental works.
AN thanks the Irish Research Council (IRC) for financial support. IR, AD and SS acknowledge additional financial support by KAUST (ACRAB project). Computational resources have been provided by Trinity Centre for High Performance Computing (TCHPC) and Irish Centre for High-End Computing (ICHEC).
- H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Nat. Phys. 5, 438 (2009).
- Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y.S. Hor, R.J. Cava, and M.Z. Hasan, Nat. Phys. 5, 398 (2009).
- D. Pesin and A.H. MacDonald, Nat. Mater. 11, 409 (2012).
- M.Z. Hasan and C.L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
- H. Beidenkopf, P. Roushan, J. Seo, L. Gorman, I. Drozdov, Y.S. Hor, R.J. Cava, and A. Yazdani, Nat. Phys. 7, 939 (2011).
- J. Wang, W. Li, P. Cheng, C. Song, T. Zhang, P. Deng, X. Chen, X. Ma, K. He, J.-F. Jia, Q.-K. Xue, and B.-F. Zhu, Phys. Rev. B 84, 235447 (2011).
- Z. Alpichshev, R.R. Biswas, A.V. Balatsky, J.G. Analytis, J.-H. Chu, I.R. Fisher, and A. Kapitulnik, Phys. Rev. Lett. 108, 206402 (2012).
- S. Kim, M. Ye, K. Kuroda, Y. Yamada, E.E. Krasovskii, E.V. Chulkov, K. Miyamoto, M. Nakatake, T. Okuda, Y. Ueda, K. Shimada, H. Namatame, M. Taniguchi, and A. Kimura, Phys. Rev. Lett. 107, 056803 (2011).
- T. Zhang, P. Cheng, X. Chen, J.-F. Jia, X. Ma, K. He, L. Wang, H. Zhang, X. Dai, Z. Fang, X. Xie, and Q.-K. Xue, Phys. Rev. Lett. 103, 266803 (2009).
- Z. Alpichshev, J.G. Analytis, J.-H. Chu, I.R. Fisher, Y.L. Chen, Z.X. Shen, A. Fang, and A. Kapitulnik, Phys. Rev. Lett. 104, 016401 (2010).
- Z. Alpichshev, J.G. Analytis, J.-H. Chu, I.R. Fisher, and A. Kapitulnik, Phys. Rev. B 84, 041104(R) (2011).
- X. Zhou, C. Fang, W.-F. Tsai, and J.P. Hu, Phys. Rev. B 80, 245317 (2009).
- R.R. Biswas and A.V. Balatsky, Phys. Rev. B 83, 075439 (2011).
- R.R. Biswas and A.V. Balatsky, Phys. Rev. B 81, 233405 (2010).
- X.-F. Wang, Y. Hu, and H. Guo, Phys. Rev. B 85, 241402(R) (2012).
- T. Nakanishi, M. Koshino, and T. Ando, Phys. Rev. B 82, 125428 (2010).
- J.M. Soler, E. Artacho, J.D. Gale, A. Garcia, J. Junquera, P. Ordejón, and D. Sánchez-Portal, J. Phys.: Condens. Matter 14, 2745 (2002).
- A. R. Rocha, V. M. Garcia-Suarez, S. Bailey, C. Lambert, J. Ferrer, and S. Sanvito, Nat. Mater. 4, 335 (2005).
- A. R. Rocha, V. M. Garcia-Suarez, S. Bailey, C. Lambert, J. Ferrer, and S. Sanvito, Phys. Rev. B 73, 085414 (2006).
- I. Rungger and S. Sanvito, Phys. Rev. B 78, 035407 (2008).
- L. Fernández-Seivane, M.A. Oliveira, S. Sanvito, and J. Ferrer, J. Phys.:Condens. Matter 18, 7999 (2006).
- J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- B. Naydenov, M. Mantega, I. Rungger, S. Sanvito, and J.J. Boland, Phys. Rev. B 84, 195321 (2011).
- M.I. Katsnelson, K.S. Novoselov, and A.K. Geim, Nat. Phys. 2, 620 (2006).
- I. A. Nechaev, M. F. Jensen, E. D. L. Rienks, V. M. Silkin, P. M. Echenique, E. V. Chulkov, and Ph. Hofmann. Phys. Rev. B 80, 113402 (2009).
- M.F. Crommie, C.P. Lutz, and D.M. Eigler, Nature 363, 524 (1993).
- I. Rungger, A. Narayan, U. Schwingenschloegl, and S. Sanvito (in preparation).
- M.L. Teague, H. Chu, F.-X. Xiu, L. He, K.-L. Wang, and N.-C. Yeh, Solid State Communications 152, 747 (2012).
- J. Seo, P. Roushan, H. Beidenkopf, Y.S. Hor, R.J. Cava, and A. Yazdani, Nature 466, 343 (2010).
- A. Narayan, I. Rungger, and S. Sanvito, Phys. Rev. B 86, 201402(R) (2012).
- The plots were generated using XCrySDen software; A. Kokalj, Comp. Mater. Sci., 28, 155 (2003).
- Y. Takane and K.-I. Imura, J. Phys. Soc. Jpn. 81, 093705 (2013).
- Y. Takane and K.-I. Imura, J. Phys. Soc. Jpn. 82, 074712 (2013).