Dirac fermions at high-index surfaces of bismuth chalcogenide topological insulator nanostructures

Binary bismuth chalcogenides Bi2Se3, Bi2Te3, and related materials are currently being extensively investigated as the reference topological insulators (TIs) due to their simple surface-state band dispersion (single Dirac cone) and relatively large bulk band gaps. Nanostructures of TIs are of particular interest as an increased surface-to-volume ratio enhances the contribution of surfaces states, meaning they are promising candidates for potential device applications. So far, the vast majority of research efforts have focused on the low-energy (0001) surfaces, which correspond to natural cleavage planes in these layered materials. However, the surfaces of low-dimensional nanostructures (nanoplatelets, nanowires, nanoribbons) inevitably involve higher-index facets. We perform a systematic ab initio investigation of the surfaces of bismuth chalcogenide TI nanostructures characterized by different crystallographic orientations, atomic structures and stoichiometric compositions. We find several stable terminations of high-index surfaces, which can be realized at different values of the chemical potential of one of the constituent elements. For the uniquely defined stoichiometric termination, the topological Dirac fermion states are shown to be strongly anisotropic with a clear dependence of Fermi velocities and spin polarization on the surface orientation. Self-doping effects and the presence of topologically trivial mid-gap states are found to characterize the non-stoichiometric surfaces. The results of our study pave the way towards experimental control of topologically protected surface states in bismuth chalcogenide nanostructures.

Binary bismuth chalcogenides Bi 2 Se 3 , Bi 2 Te 3 , and related materials are currently being extensively investigated as the reference topological insulators (TIs) due to their simple surface-state band dispersion (single Dirac cone) and relatively large bulk band gaps. Nanostructures of TIs are of particular interest as an increased surface-to-volume ratio enhances the contribution of surfaces states, meaning they are promising candidates for potential device applications. So far, the vast majority of research efforts have focused on the low-energy (0001) surfaces, which correspond to natural cleavage planes in these layered materials. However, the surfaces of low-dimensional nanostructures (nanoplatelets, nanowires, nanoribbons) inevitably involve higher-index facets. We perform a systematic ab initio investigation of the surfaces of bismuth chalcogenide TI nanostructures characterized by different crystallographic orientations, atomic structures and stoichiometric compositions. We find several stable terminations of high-index surfaces, which can be realized at different values of the chemical potential of one of the constituent elements. For the uniquely defined stoichiometric termination, the topological Dirac fermion states are shown to be strongly anisotropic with a clear dependence of Fermi velocities and spin polarization on the surface orientation. Self-doping effects and the presence of topologically trivial mid-gap states are found to characterize the non-stoichiometric surfaces. The results of our study pave the way towards experimental control of topologically protected surface states in bismuth chalcogenide nanostructures.
Topological insulators (TIs) are a recently discovered class of materials that realize a non-trivial band structure topology induced by strong spin-orbit interactions [1][2][3] . These novel materials possess topologically protected metallic surface states characterized by robust spin-momentum locking. Initial theoretical predictions [4][5][6][7] and subsequent experimental realizations of bulk TIs 8 have generated an enormous degree of interest in exploring this exotic electronic phase. Significant attention has focused on the so-called "second generation" of bulk TIslayered binary bismuth chalcogenides Bi 2 X 3 (X = Se or Te) 9,10 . The simple dispersion of their surface-state bands, represented by a single Dirac cone, and relatively large band gaps (~0.3 eV for Bi 2 Se 3 ), indicating room temperature stability of the topological phase, have made these materials the "reference" TIs. Despite significant progress in probing surface states via surface sensitive techniques such as angle-resolved photoemission spectroscopy 9,11,12 (ARPES) and scanning tunnelling microscopy/spectroscopy [13][14][15] (STM/STS), electron transport measurements have proven to be more challenging. This is primarily due to surface-state charge carriers being masked by bulk states as a result of intrinsic impurities, such as vacancy and antisite defects 16,17 , or the ambient environment 18,19 .
An efficient strategy towards mitigating bulk contributions is to form low-dimensional nanostructures with an increased surface-to-volume ratio [20][21][22][23][24] . Synthesized bismuth chalcogenide nanostructures display a diverse array of morphologies such as hexagonal and triangular platelets, nanowires and nanoribbons, including both linear and irregularly shaped nanoribbons (for example see ref. 20). In addition to giving access to transport measurements of topologically protected surface states [25][26][27][28] , low-dimensional TI nanostructures also provide the necessary configuration for observing a variety of physical phenomena, for example the Aharonov-Bohm interference [29][30][31] . Finally, one can anticipate novel physical phenomena involving topologically protected charge carriers subjected to reduced dimensions, such as the realization of Majorana fermion quasiparticles 32 .
Scientific RepoRts | 6:20220 | DOI: 10.1038/srep20220 Bismuth chalcogenides are layered materials composed of covalently bonded sheets termed quintuple layers (QL) (Fig. 1a-c). These QL building blocks are held together by weak van der Waals (vdW) interactions, thus providing natural cleavage planes and giving rise to two equivalent orientations of low-energy surfaces, (0001) and (0001). The vast majority of investigations of the surface states of bismuth chalcogenide TIs thus far have focused on the low-energy (0001) surfaces. However, it has to be appreciated that any bismuth chalcogenide nanostructure of dimensionality lower than 2 inevitably exhibits surfaces with orientations other than the two abovementioned. The properties of topologically protected charge carriers at such surfaces would depend on their crystallographic orientation, atomic structure and chemical composition. Furthermore, the properties of a TI nanostructure as a whole are defined by the relative presence of different facets on the nanostructure's surface, which is in turn determined by their respective surface energies.
Fully realizing the potential of bismuth chalcogenide nanostructures for exploring the fundamental physics and device applications of TIs depends on a clear understanding of how structure and morphology of the surfaces impact upon the electronic properties of those nanostructures. So far, only one example of a high-index surface of bismuth chalcogenides has been addressed theoretically 33 and experimentally 34,35 . In this paper, we report a systematic first-principles investigation of high-index surfaces of Bi 2 Se 3 and Bi 2 Te 3 . We first determine the structure and chemical composition of stable surfaces. Subsequently, the electronic properties of corresponding topological surface-state charge carriers and their dependence on surface orientation and local chemical composition are investigated.

Results
Stable quintuple-layer terminations and surface energies. Any investigation of the electronic structure of a material's surface requires detailed knowledge of the atomic structure of stable configurations of that surface. In layered materials, such as Bi 2 Se 3 and Bi 2 Te 3 , the surface structure is defined by the termination of individual layers, and by their orientation with respect to the surface plane, which we here define in terms of the angle θ, as shown in Fig. 1c. In this definition, θ = 0° corresponds to the degenerate case of the low-energy (0001) surface, for which no QL termination takes place. Furthermore, given that Bi 2 Se 3 and Bi 2 Te 3 are layered vdW materials, the surface energy E (expressed in eV/Å 2 ) can subsequently, to a good degree of accuracy, be related to the termination energy of individual quintuple layers G (expressed in eV/Å) as where c is the lattice constant of a hexagonal unit cell of the bismuth chalcogenide material, and hence c/3 is the QL thickness (9.546 Å for Bi 2 Se 3 and 10.162 Å for Bi 2 Te 3 36 ). The approximation in Equation 1 stems from the omission of a term related to the loss of vdW energy upon stacking QLs to form the 2D surface. This contribution is assumed to be negligible as the relative magnitude of vdW interactions is small and the leading term to the surface energy is the QL termination energy. Below, as we will be referring to the hexagonal unit cell of bismuth chalcogenides, correspondingly, crystallographic planes and lattice directions will be given in the four-index Miller-Bravais notation. One has to appreciate the fact that QL terminations can have different local stoichiometries, either Bi-rich or Se(Te)-rich, and that their relative stabilities depend on the chemical potential of one of the constituent elements, μ Bi or equivalently μ Se (μ Te ). The chemical potential reflects experimental conditions under which the nanostructure is grown, and, importantly, its varying values may result in different stable surface structures. Hence, it is possible to tailor the structure of high-index surfaces by changing the experimental (c) Two-dimensional slab model illustrating the structure of one example of a high-index surface (left facet). The structure of high-index surfaces of layered bismuth chalcogenides is defined by stacking angle θ and the QL termination. The example shown corresponds to θ = 57.7° and the stoichiometric QL termination (configuration I).
conditions, which consequently translates into control over the electronic properties of topologically protected states hosted by these surfaces.
The energies of QL terminations have been systematically investigated by means of density functional theory (DFT) calculations carried out on a large number of single-QL models in a nanoribbon configuration (see Methods). Relative stabilities of different QL terminations were determined by comparing the QL termination free energies G(μ Bi ), which are calculated per unit length as a function of the chemical potential of Bi, µ Bi , via the expression Bi Bi with G(μ Bi = 0) being the QL termination free energy at a reference chemical potential where X corresponds to Se (Te). Here, G(μ Bi = 0) is calculated as the formation energy of a QL termination defined as the difference in total energies, per unit length, of a nanoribbon model (E model ), isolated two-dimensional QL per Bi 2 Se 3 (Bi 2 Te 3 ) unit ( ) E Bi X 2 3 , and Bi atom in its bulk elemental crystal (E Bi ). The reference chemical potential μ Bi = 0 eV below corresponds to bulk elemental bismuth. N Bi refers to the number of excess Bi atoms in the model relative to a stoichiometric system (i.e. with Bi : X = 2 : 3), while N Bi X 2 3 is the number of stoichiometric Bi 2 X 3 units. L is the periodicity of the nanoribbon model, equal to the lattice constant a of hexagonal unit cell (4.178 Å for Bi 2 Se 3 and 4.403 Å for Bi 2 Te 3 , which correspond to values determined here from first-principles) in the case of the QL edge oriented along the [2110] direction, and to = L a 3 in the case of the QL edge oriented along the [1100] direction. The [2110] and [1100] lattice directions are labelled in reference to the hexagonal unit cell in Fig. 1b.
The calculated QL termination free energies G(μ Bi ) in Bi 2 Se 3 and Bi 2 Te 3 are presented in Fig. 2a,b, respectively. The dependence of G(μ Bi ) on the chemical potential μ Bi is linear, with the slope reflecting the local deviation from the stoichiometric ratio Bi : X = 2 : 3 at the QL edge. Thicker coloured lines indicate stable terminations, i.e. those having the lowest value of G(μ Bi ) within a certain range of μ Bi . For Bi 2 Se 3 , we found 7 such configurations: 1 stoichiometric (labelled I in Fig. 2a), 3 Se-rich (II Se , III Se and IV Se ) and 3 Bi-rich (II Bi , III Bi and IV Bi ). Their atomic structure configurations are schematically shown in Fig. 2c. The uniquely defined stoichiometric termination, characterized by the constant G = 0.169 eV/Å, is stable in a particularly wide range of μ Bi values, and it displays only very minor structural relaxation (see Fig. S1 of the Supplementary Information). Several distinct non-stoichiometric terminations are possible because of the varying deviations from ideal stoichiometry indicated in Fig. 2c. We note that the stable terminations tend to be oriented along the [2110] direction, with only the II Se configuration being oriented along the [1100] direction. This implies that high-index surfaces formed by the [2110] QL terminations are preferred from the thermodynamic point of view. On the other hand, an observation of a surface formed by the [1100] QL termination unambiguously points to the structure and stoichiometry of the II Se configuration. Analysis of the relaxed atomic structures (Fig. S1 of the Supplementary Information) suggests that the non-stoichiometric configurations with the largest deviations from the stoichiometric ratio, IV Bi and IV Se , can be interpreted as the stoichiometric edge configurations I capped by elemental Bi and Se, respectively. This reflects the expected tendency towards the onset of phase segregation into Bi 2 Se 3 and elemental Bi and Se, respectively, at the extreme values of chemical potential μ Bi . For Bi 2 Te 3 , we predict very similar structures of stable configurations and their associated energies (Fig. 2b). For instance, the stable stoichiometric configuration I with G = 0.161 eV/Å also exists in a broad range of μ Bi values. However, we find that the QL termination analogous to III Se is not present on its respective phase diagram.
The knowledge of surface energies E for different crystallographic orientations of the surface allows one to deduce the equilibrium crystal shape by means of the Wulff construction 37 . In other words, it enables the contributions of different facets to the overall surface of the nanostructure to be determined. The surface energies of high-index facets can be deduced from the QL termination energies discussed above, except for the (0001) orientation, which is of the vdW origin. We estimate the (0001) surface energy as E vdW = 0.01 eV/Å 2 , using the recently suggested universal interlayer binding energy of layered materials 38 , and a first-principles study of the role of vdW interactions in the binding energy of Bi 2 Se 3 and Bi 2 Te 3 39 . This value can be compared with energies of high-index surfaces for θ = 90° by plotting G equiv = E vdW c/3 (dashed line in Fig. 2a,b). In the range of chemical potentials where the most stable QL termination is the stoichiometric configuration I, the low-energy (0001) surface will be dominant. Assuming thermodynamic equilibrium conditions of nanoparticle growth, this will lead to the nanoplatelet or nanoribbon morphology with top and bottom facets being the (0001) and (0001) low-energy surfaces, while the side surfaces would have the stoichiometric composition with the QL terminations corresponding to configuration I (Fig. 2d). Such nanoplatelets of bismuth chalcogenides have been extensively documented in experimental publications 21,22,28 . For chemical potential values corresponding to Bi-rich or Se(Te)-rich conditions, the energies of stable non-stoichiometric surfaces appear to be lower than those of the (0001) and configuration I stoichiometric surfaces. This implies the nanowire morphology for nanostructures 20,31 grown under such conditions, with the surface dominated by the non-stoichiometric compositions and oriented along the (0001) direction, for example the hexagonal cross-section nanowires shown schematically on the right and left hand side of Fig. 2d. From a general perspective, this points to a possibility of tailoring the morphology, composition and properties of bismuth chalcogenide nanostructures by controlling the growth conditions. Scientific RepoRts | 6:20220 | DOI: 10.1038/srep20220 Topologically protected states at stoichiometric high-index surfaces. We now turn our discussion to the electronic structure of high-index surfaces of bismuth chalcogenide TIs. As discussed above, the atomic structure of such surfaces is determined by the crystallographic orientation of the surface and the termination of individual QLs. Initially, we focus on the first degree of freedom assuming the stoichiometric termination I, which is stable in wide ranges of chemical potentials in both materials (Fig. 2a-c). Calculations have been performed on two-dimensional slab models in which the interior part reproduces bulk crystal structure (see Methods). One example is shown in Fig. 1c, and corresponds to θ = 57.7° and θ = 58.0° for Bi 2 Se 3 and Bi 2 Te 3 , respectively, where θ defines the QL stacking angle of a high-index surface. In this particular surface configuration, atomic planes along which the QLs are terminated almost coincide with the surface place. Figure 3a   orientation-dependent band dispersion of surface-localized Dirac fermion states, which agrees well with the results of previous calculations performed on this particular surface 33 -the only high-index surface of Bi 2 Se 3 addressed theoretically. More importantly, this is also the only case of a high-index surface that was investigated experimentally using samples of Bi 2 Se 3 epitaxially grown on the InP(001) substrate 34,35 . The degree of Dirac fermion anisotropy, quantified as the ratio of Fermi velocities for the momenta along x′ direction (i.e. along the QL planes, see Fig. 1c) and along y′ direction (perpendicular to x′ ), υ υ / ′ ′ x y F F = 2.07, agrees well with υ υ / ′ ′ x y F F ~ 2 observed in the ARPES measurements of ref. 34. Interestingly, the same surface configuration of Bi 2 Te 3 shows an even more pronounced anisotropy of the Dirac fermion surface states with practically flat bands for the momenta along the y′ direction (Fig. 3c).
Band structures have been calculated for other crystallographic orientations of the surface within the range 25° < θ < 155°. Three representative surface configurations are shown in Fig. 4a. For all studied surface orientations the band structures feature an anisotropic Dirac cone at the Γ point, but the degree of anisotropy varies across the investigated range of θ. Figure 4b shows Fermi velocities computed above the Dirac point energy and for momenta along x′ and y′ as a function of surface orientation θ, for both Bi 2 Se 3 and Bi 2 Te 3 . One can notice the following systematic trends in the above-mentioned quantities. Firstly, Fermi velocities of the topological-state bands of high-index surfaces are generally lower than that for the (0001) surface (4.8 × 10 5 m/s and 4.6 × 10 5 m/s for Bi 2 Se 3 and Bi 2 Te 3 , respectively, obtained from our calculations). Secondly, the largest anisotropies are achieved around θ = π 2 , which corresponds to QL planes oriented perpendicular to the surface. In the case of Bi 2 Se 3 , the value of υ υ / ′ ′ x y F F does not exceed 2.5, while for Bi 2 Te 3 the surface-state band is practically dispersionless along y′ in a broad range of values of θ resulting in much larger anisotropies.
Spin texture is another important property of topologically protected surface states. This property can be analysed by investigating the momentum dependence of the expectation values of spin operators ). Figure 4c shows the magnitude of the spin-polarization vector P of the surface-state charge carriers above the Dirac point energy as a function of surface orientation θ for Bi 2 Se 3 and Bi 2 Te 3 . The magnitude of spin-polarization P is always lower along the x′ direction, with the largest difference achieved for surface orientations around θ = π 2 , similar to the case of the Fermi velocities. In contrast, the average magnitude of spin-polarization is about the same as for the (0001) surface. This can be compared to the model theory results of refs 44,45, which predict that at θ = π 2 the spin texture collapses to a single dimension along the x′ direction, leading to a complete suppression of spin polarization for momenta oriented along this direction. In our first-principles calculations, however, we do not observe such a vanishing spin polarization, nonetheless its magnitude for momenta along the x′ direction achieves its minimum down to ~0.4 for both chalcogenide materials.
Electronic structure of non-stoichiometric high-index surfaces. The models of non-stoichiometric surfaces are found to reveal generally more complex band structures. In particular, in addition to the topologically protected crossing within the band gap, we observe the presence of topologically trivial mid-gap states and self-doping, i.e. charge transfer between the surface and bulk-like states that leads to the shift of the Fermi level into the valence and conduction bands. Both effects can be considered as detrimental to the observation of topological surfaces states in TI nanostructures. Figure 5 shows two representative examples of band structures of non-stoichiometric surface models. In the case of a Se-rich surface of Bi 2 Se 3 (III Se termination, Fig. 5a) we observe the presence of topologically trivial states with high band dispersion crossing the band gap and strongly hybridizing with the topologically protected states. The topologically trivial states are distinguished from the topologically  protected, the presence of which is guaranteed by the strong topological nature of Bi 2 Se 3 and Bi 2 Te 3 , by means of the bulk-boundary correspondence. The bands corresponding to the trivial states cross some energy within the band gap an even number of times, while an odd number of crossings can be observed for the topologically protected states. In the case of a Bi-rich surface of Bi 2 Se 3 (II Bi termination, Fig. 5b) the Fermi level appears shifted into the conduction band, i.e. n-type doping of the bulk-like states is observed. At the same time, the Dirac point of the topological surface-state band at the Γ point appears to be immersed in the valence band. The change in the position of the Dirac point at Γ pulls another topologically protected crossing at the time-reversal invariant momentum point Y into the bulk band gap. A similar reorganization of the band dispersion of the topological surface states upon changing the surface structure has been predicted for a model bulk topological insulator 46 . An overview of the general electronic structure features of non-stoichiometric surfaces of Bi 2 Se 3 and Bi 2 Te 3 are given in Table 1, while the corresponding band structures are reproduced in Figs S2 and S3 of the Supplementary Information. In general, we find that n-type doping is a characteristic feature of all Bi-rich surfaces, while the presence of trivial mid-gap states is observed for both chalcogen and Bi-rich surfaces in the extreme limits of chemical potential μ Bi .

Discussion
In summary, our work addresses the atomic structure and electronic properties of high-index surfaces in nanostructures of bismuth chalcogenide topological insulators, Bi 2 Se 3 and Bi 2 Te 3 , through first-principles calculations. We predict that several possible quintuple-layer terminations of different stoichiometric compositions can be realized, depending on experimental conditions. Both the stoichiometry of the surface and its crystallographic orientation significantly affect the electronic properties of topologically protected surface states, particularly the anisotropy of their Dirac fermion band dispersion and the degree of spin polarization. Moreover, these properties are shown to display clear dependence on the surface configuration. Through this understanding, one can gain a greater degree of control over the properties of nanostructures of topological insulators aiming at prospective technological applications of these novel materials.

Methods
First principles calculations were performed within the DFT framework, employing the generalized gradient approximation (GGA) to the exchange-correlation functional 47 . Spin-orbit effects were treated self-consistently using fully relativistic norm-conserving pseudopotentials 48 within the two-component wavefunction formalism. A plane-wave kinetic energy cutoff of 35 Ry has been employed for the wave functions. Calculations were implemented through the PWSCF plane-wave pseudopotential code of the Quantum-ESPRESSO distribution 49 .  First-principles band structures of slab models of (a) Se-rich (III Se QL termination) and (b) Bi-rich (II Bi QL termination) high-index surfaces of Bi 2 Se 3 at θ = 57.7°. Points A and Y correspond to the Brillouin zone boundaries along directions defined by reciprocal lattice vectors associated with the real-space unit vectors of the surface. The size of red symbols reflects the magnitude of the inverse participation ratio (IPR).
Scientific RepoRts | 6:20220 | DOI: 10.1038/srep20220 Quintuple-layer termination energies at reference chemical potential G(μ Bi = 0) were calculated using single-QL nanoribbon models of 5.4 nm and 5.7 nm width for Bi 2 Se 3 and Bi 2 Te 3 , respectively. In total, over 60 different termination structures were considered, shown as grey lines in Fig. 2a,b. The electronic structure of high-index surfaces was investigated using two-dimensional slab models of the same thickness as the width of nanoribbon models. The lattice constants of slab configurations were calculated assuming bulk lattice constants as well as the stacking order of bulk crystals.