Topological surface states on Bi(111) based on empirical tight-binding calculations

The topological order of single-crystal Bi and its surface states on the (111) surface are studied in detail based on empirical tight-binding (TB) calculations. New TB parameters are presented that are used to calculate the surface states of semi-infinite single-crystal Bi(111), which agree with the experimental angle-resolved photoelectron spectroscopy results. The influence of the crystal lattice distortion is surveyed and a topological phase transition is found that is driven by in-plane expansion. In contrast with the semi-infinite system, the surface-state dispersions on finite-thickness slabs are non-trivial irrespective of the bulk topological order. The role of the interaction between the top and bottom surfaces in the slab is systematically studied, and it is revealed that a very thick slab is required to properly obtain the bulk topological order of Bi from the (111) surface state: above 150 biatomic layers in this case.


I. INTRODUCTION
Topological materials classified by unconventional parity eigenvalues of three-dimensional (bulk) bands is one of the topics of high interest in solid state physics in this decade.The insulators classified in the non-trivial (topological) group are called topological insulators (TIs), and hold metallic and spin-polarized surface states that continuously disperse between the bulk valence band maximum and the conduction band minimum [1][2][3].Because these topological surface states are spin-polarized and robust against any perturbations that do not change the parities of the entire bulk band structure, these states are regarded as a promising element for future spintronic devices [4].
In the early stages of TI research, first-principles calculations based on density functional theory (DFT) achieved great success in predicting the topology of many materials and in proposing new TI candidates [5][6][7][8][9].Most of these predictions were soon proven experimentally and there was excellent agreement between the predicted topological surface states and the observations [10][11][12].However, despite these great successes, there remains an open question on the topological order of the very simple material of single-crystal Bi.
Bismuth is widely used as a component of TIs, such as Bi 2 Se 3 [5,10] and TlBiSe 2 [7][8][9]12], because it is the heaviest non-radioactive element that also possesses a strong spin-orbit interaction (SOI).This is critical because the SOI plays an important role to realize the unconventional parity eigenvalues.The topological order of single-crystal Bi has also been extensively studied together with its relative alloy Bi 1−x Sb x , which is the first material experimentally detected as three-dimensional TI [13,14].According to the DFT [1,13] as well as empirical tight-binding (TB) [15,16] calculations, the topological order of Bi is trivial, and alloying with Sb causes the topological phase transition to a TI to occur at x ≈ 0.04 .However, unlike most of the other cases, this prediction is not consistent with the experimental results.
Figures 1(a) and 1(b) show the experimental surface-band dispersion of Bi(111) obtained by angle-resolved photoelectron spectroscopy (ARPES) [17].In these dispersions, the two spin-split surface bands S1 and S2 merge into the same projected bulk valence bands (BVBs) at Γ.However, at M , these spin-split surface bands separately merge into two different projected bulk bands, with S1 merging into the bulk conduction band (BCB), while S2 merges into the BVB.The overall dispersion is also schematically depicted in Fig. 1 on this surface-state dispersion, S1 continuously connects the projected BVB and BCB between two time-reversal-invariant momenta (TRIM), which is the behavior that is expected for topologically non-trivial materials.Indeed, a recent ARPES report on Bi (1−x) Sb x (x ∼ 0.1) has demonstrated an almost identical surface-state dispersion [18], with the sole difference being that the projected BCB is above the Fermi level at M in Bi (1−x) Sb x .It should be noted that the topological classification of the band structure is not only valid for semiconductors, but also for semimetals whose finite projected bulk band gap opens in any k in the surface plane.Interestingly, the surface-band dispersions observed by ARPES are qualitatively the same in various experiments on single-crystal Bi(111) [19] as well as thin films possessing a few tens of biatomic layers (BL) grown on Si(111) [18,20].Specifically, all the observations show that S1 merges into the BCB while S2 merges into the BVB at M .The surface-state dispersions simulated by theoretical calculations depend upon the computational methods, however.The DFT calculations based on the local density approximation (LDA) with SOI, using slab geometry to mimic the crystal surface, obtain the surface bands depicted in Fig. 1(e) [20,21].In this case, both S1 and S2 merge into the BVB at M , and hence the surface band dispersion is trivial.This is consistent with the theoretical prediction mentioned above, but disagrees with the ARPES experiments.The other major method used to calculate the surface state is the so-called transfer-matrix (TM) method with the empirical TB model for bulk electronic states [16,22].Based on this method, however, the calculated surface state exhibits an additional crossing between Γ and M, as shown in Fig. 1(f).This unexpected crossing can be understood as the influence arising from an incorrect mirror Chern number via the empirical TB parameters [16].Further, another TM calculation based on a different set of TB parameters [23] has resulted in surface bands without this unrealistic surface-state crossing, and is shown in Fig. 1(g).In this case, the Dirac point lies in the projected bulk band gap at M .The result in Fig. 1(g) is also a topologically trivial surface-band dispersion because no surface band continuously connects the BVB and BCB.

(d) . Based
The main reason for the difficulty in calculating the electronic structure of single-crystal Bi is that its bandgap is very small.The bandgap for single-crystal Bi is ∼15 meV at L in the bulk Brillouin zone, corresponding to M on the surface Brillouin zone as shown in Fig. 1(c).One of the well-known weaknesses of DFT is its inability to estimate the accurate size of the bandgap.Actually, LDA overestimates the size of the bandgap at L [21].A recent study based on the quasiparticle self-consistent GW method including SOI has improved the size of the bandgap [24].Even using such state-of-the-art computational methods, however, the topological order of single-crystal Bi is still calculated to be trivial, which disagrees with the experimental results.In the TB calculation, the size of the bandgap agrees with the experiments because the TB parameters are empirically tuned to reproduce these experimental results, although the topological order was also calculated to be trivial.
The tiny bulk bandgap at L leads to "flagile" topological phase of single-crystal Bi, because various perturbation from strain, ultrathin film thickness and so on can invert the energetical order of the bulk bands at L. Actually, a DFT calculation taking the intersurface interaction in the slab into account showed non-trivial surface-band dispersion as shown in Fig. 1(d) [20].Recently, a TB calculation using the slab model also showed the surface states which disperses from BVB at Γ to BCB at M continuously [25].Since the bulk topological order based on DFT and TB with known parameter set is trivial, this result suggests the topological phase transition driven by ultrathin film thickness of Bi.However, the magnitude of such finite-size effect, in other words, how thick the slab should be in order to calculate the surface states of Bi obeying the bulk topological order, has not been studied yet.Structural strain is also claimed as a source of topological phase transition.It is claimed that the in-plane structural strain in Bi(111) ultrathin film causes the topological phase transition from trivial to non-trivial phase [23,26].A recent theoretical study supports this results [24].However, there is still a discontinuity to the bulk case without strain: based on the ARPES experiments, Bi without strain should be topologically non-trivial and hence it is not clear where the topological phase transition takes place with lattice distortion.
In this work, we present the new set of TB parameters for calculating the surface states of single-crystal Bi(111), which agrees with the experimental ARPES results.Two sets of parameters were obtained that generate topologically trivial and non-trivial surface states for Bi, where neither exhibit any artificial crossing such as that generated previously (Fig. 1(e)) [15,16].Based on the new TB parameters, the surface-state calculations were performed by means of both the TM method and slab geometry.In both cases, the calculated surface states agree well with previous experimental results, except for in the proximity of M.
Around M , the surface-state dispersion changes depending upon the topological order of the bulk bands, where only the topologically non-trivial case agrees with the experiments.The influence of crystal lattice distortion was surveyed and a topological phase transition was found that was driven by in-plane expansion for the non-trivial bulk bands, which was opposite to the distortion suggested in a previous report [23].In contrast with the semi-infinite crystal, the slab calculation generated non-trivial surface-state dispersions irrespective of the bulk topological order, as has been the case for previous DFT calculations using slab geometries.The role of the interaction between the top and bottom surfaces in the slab was systematically studied, and it was found that a very thick slab is required to obtain the bulk topological order of Bi from the (111) surface states properly.Specifically, the slab must be greater than 150 BL in our model.

A. Tight-binding parameters for bulk states
The main framework used to calculate the bulk electronic structure was the same as that in Ref. 15 , and is briefly explained herein.Single-crystal Bi has an A17 rhombohedral lattice, but is also characterized by hexagonal lattice parameters, a and c, together with an additional parameter µ (see Fig. 2).The primitive (rhombohedral) unit cell contains two atoms and the relative position of the second atom is (0, 0, 2µc), where µ = 0.2341.The hopping parameters of the sp 3 orbitals between first-, second-and third-nearest-neighbor atoms were taken, and the SOI was included via the spin-orbit coupling parameter λ = 1.5 eV.The resulting (16×16) matrix is shown in the appendix of Ref. 15 .
Next, the TB parameters were modified so that the calculated surface-state dispersion between Γ and M could be reproduced without any artificial crossing such as was obtained in Ref. 16 .In addition, the new parameters were tuned to maintain the energies of the electron/hole pockets and the size of the bandgap at L at nearly the same value as the original parameter, which agrees with the experimental values.II, the major difference between TBP-1 and TBP-2 is the difference of topological order: trivial for TBP-1 but non-trivial for TBP-2.

B. Transfer-matrix (TM) method
The TM method is used to calculate the surface electronic states on a semi-infinite crystal from the bulk Hamiltonian and the transfer matrix T (k , E), as proposed in Ref. 22 .In this TM, k is the in-plane wavevector and E is the binding energy .The procedure reported in a previous paper [16] was followed, as described below, but the bulk TB parameters were and the Z 2 topological invariants (ν 0 ; ν 1 ν 2 ν 3 ) calculated with the tight-binding parameters shown in Table I.
Because the Bi crystal can be regarded as a stack of BLs, as depicted in Fig. 2(b), the bulk electronic states of a semi-infinite Bi crystal can be written as where φ n,a is a basis of the states in the BL plane localized on the a = 1, 2 monolayer of the nth BL.In addition, each φ n,a has eight components associated with the eight atomic orbitals.The transfer matrix, T (k , E), is given by Equations

C. Finite slab calculation
Slab geometry is widely used for DFT calculations of surface electronic structures, such as used in Refs.20, 21, 23, 25 .This model can mimic the surface without breaking the three-dimensional periodicity and can therefore be easily applied toward various surface systems.However, sometimes this model generates artificial states owing to the interaction between the top and bottom surfaces.In this work, we followed the method reported in the recent paper [25], as is described below.Varying from the previous work, however, we used a new set of bulk TB parameters and assumed no surface hopping term.
The total Hamiltonian of the slab H is represented as 12 H 11 H (1) 12 where H 11 is the hopping term for inter-monatomic-layer hopping (i.e., in the same monatomic layer) and H (1) 12 (H 21 ) are the intra-BL (inter-BL) hopping terms.The H 12 , H 21 , and H 11 terms are given in the bulk TB Hamiltonian [15] as the first-, second-and thirdnearest-neighbor hopping terms, respectively.The size of the matrix is therefore 16n×16n, where n is the number of the BLs in the slab.
In the following, the obtained states were plotted with Gaussian broadening (F W HM = 5 meV) to mimic the ARPES intensity plots produced with a typical instrumental energy resolution.

III. RESULTS AND DISCUSSION
A. Surface states generated by TM method The same as (a) and (b), respectively, but calculated using the new TB parameters generating topologically trivial bulk bands (TBP-1 in Table I).(e,f) The same as (a) and (b), respectively, but calculated using the new TB parameters generating topologically non-trivial bulk bands (TBP-2 in Table I).
band dispersion qualitatively agrees with the ARPES results except for the k region around M , as shown in Figs.3(c) and 3(d).The two branches of the surface bands degenerate with each other at M , and this surface-state dispersion is therefore topologically trivial, as is expected from the topological order of bulk bands.This surface-state dispersion agrees with that reported in a previous paper [23].Around Γ, the TB parameters generating nontrivial bulk bands (TBP-2 in Table I) does not significantly alter the surface-state dispersion from the trivial dispersion, as shown in Fig. 3(e).Around M ,however, the surface-state dispersion is different from those calculated with the other TB parameters.Specifically, the upper branch merges into the BCB while the lower merges into the BVB, showing a good agreement with the experimental results [17][18][19][20] (see Fig. 3(f)).
These good agreements of surface-state dispersions with the previous experimental and theoretical results except for the proximity of M implies that the TB calculation can no longer be the evidence of topological order of single-crystal Bi.In order to judge such balancing two possibilities, one requires the experimental data.Based on the ARPES results, it should be topologically non-trivial, while ARPES does not provide any direct information about the parity eigenvalues of bulk bands at L. The other, bulk-sensitive and accurate experimental method would be helpful to make a firm conclution on this controversial issue, topological order of single-crystal Bi.

B. Topological phase transition driven by lattice distortion
In order to examine the topological phase transition driven by structural distortion based on Refs.23, 24, 26, we surveyed the topological order and surface-state dispersion around M using the TM calculation and the new TB parameters.
Figure 4(a) shows the bulk band evolution at L with an in-plane lattice distortion inserted for calculations using the two TB parameter sets in Table I, TBP-1 and TBP-2.Irrespective of the topological order generated with zero lattice distortion, a topological phase transition occurs in both cases.The only difference observed is whether the topological phase transition occurs when one uses a lattice strain (TBP-1, trivial at zero lattice distortion) or a lattice expansion (TBP-2, non-trivial at zero lattice distortion).Figures 4(b-e) show the surfacestate dispersions around M obtained using the TM calculation as Fig. 3.Note that the M position of each plot changes according to the in-plane lattice constant.The plots in Figs.

4(b) and 4(d)
show that the surface-state dispersion is topologically non-trivial with lattice strain (−5 %) for both TB parameter sets.In addition, with an in-plane lattice expansion (+5 %), the surface states for both TB parameter sets exhibit a trivial dispersion, as seen in Figs.4(c) and 4(e), that form a Kramers-degeneracy point at M in the projected bulk  I.
bandgap.It should be noted that the non-trivial surface-state dispersion observed by ARPES experiments with the presence of in-plane lattice strain [23,26] exhibits no conflict with the calculated results from both of the new TB parameters, TBP-1 and TBP-2.
Based on these results, we propose that to trace the bulk bandgap of single-crystal Bi with in-plane lattice distortion in to conclude the order of Bi.If the bandgap closes with the in-plane lattice expantion (strain), it would be the smoking-gun evidence of non-trivial (trivial) topological order.

C. Surface states on a finite slab
Figure 5 shows the electronic structure calculated with the slab geometry, where the slab thickness is 20 BL and the dashed lines represent the edge of the projected bulk bands.The surface-state bands dispersing in the projected bulk bandgap are obtained together with the discrete quantum well (QW) states corresponding to the bulk bands in the projected bulk-band region.Around Γ, the surface-state dispersions are almost identical to those calculated by the TM method (see Figs. 3(c) and 3(e)).Note that no additional surface hopping terms are needed to obtain these surface bands, which is in contrast with a previous study [25].The fact that no additional hopping terms are needed is possibly owing to the different TB parameter set used in this study; however, the surface-state dispersions are quite different around M. Even when using the TB parameter that generates a trivial topology of bulk bands (TBP-1 in Table I), the surface-state dispersion obtained by the slab geometry suggests a non-trivial topological order wherein the upper branch merges into the BCB while the lower merges into the BVB.Such behavior is the same as that reported in previous studies [20,25], wherein it has been explained as the influence of the interaction between the top and bottom surfaces owing to the finite slab thickness.Even when using the trivial TB parameter set (TBP-1 in Table I), the inter-surface interaction would re-invert the bulk bands at L, where the bulk bandgap is the smallest, and thus cause the topological phase transition.Similar topological phase transitions driven by the film thickness have been observed in the QWs of HgTe [27,28].
To explore the finite-thickness effect in more detail, we calculated the electronic states around M at different thicknesses.Figures 6(a) and 6(b) plot the electronic structure around M in a 200 BL slab calculated with the TB parameters TBP-1 and TBP-2, respectively, which are given in Talbe I.The number of the QW states in these plots are much larger than those of the 20-BL slab (cf. Figure 5), reflecting the increased slab thickness.In the proximity of M , however, one can find the difference between Figs. 6(a) and 6(b).The two surface-state branches dispersing out of the projected bulk bands in Fig. 6(b) enter the projected bulk bands and become QW states.However, in Fig. 6(a), these branches do not enter the projected bulk bands but remain in the bulk bandgap at M , as is the case for the semi-infinite TM calculation (cf.Fig. 3(d)).Figure 6(c) plots the energy evolution of these two surface states at M , which are connected to the surface states away from M , together with the energy positions of the projected BVB and BCB.As shown in Fig. 6(c), the surface states calculated with TBP-1 (trivial bulk bands) disperse out of the projected bulk bands at a slab thickness greater than ∼150 BL.In contrast, with TBP-2 (topologically non-trivial bulk bands), the surface states never appear out of the projected bulk bands at M .This result suggests that a much thicker slab is required to accurately measure the topological order of single-crystal Bi than has been used in previous studies [23,26].With a slab that is insufficiently thick, the surface-state dispersion always exhibits the topologically non-trivial behavior irrespective of the topological order of the bulk Bi.
Recently, two groups have reported ARPES experimental results of the Bi(111) surface states using a thickness above 100 BL [29,30].In both cases, the surface-state bands disperse as that seen in Fig. 6(b), and hence these experimental results indicate the nontrivial topological order of Bi, which is in good agreement with the bulk single-crystal case [17].

IV. SUMMARY
In summary, new TB parameters are presented with which to calculate the surface states of single-crystal Bi(111) that agree with the experimental ARPES results.Two sets of TB parameters were obtained that make Bi topologically trivial or non-trivial, wherein neither parameter set exhibits any artificial crossings such as those generated by previous TB parameters (Fig. 1(e)) [15,16].Based on the new TB parameters, surface-state calculations were performed using both the transfer-matrix method and slab geometry.In both cases, the calculated surface states agree well with the previous experimental results, except for in the proximity of M .Around M , the surface-state dispersion changes depending upon the topological order of the bulk bands, wherein only the topologically non-trivial case agrees with the experiments.We surveyed the influence of crystal lattice distortions and found a topological phase transition driven by in-plane expansion for the non-trivial bulk bands, which is opposite to the distortion found in a previous report [23].In contrast with a semi-infinite crystal, the slab calculation generated non-trivial surface-state dispersions irrespective of the bulk topological order, as is the case suggested in previous DFT calculations using slab geometries.The role played by the interaction between the top and bottom surfaces in the slab was systematically studied and it was revealed that a very thick slab (i.e., greater than 150 BL in our case) is required to accurately obtain the bulk topological order of Bi from the (111) surface states.These detailed calculations of the electronic structure and topological order of the simple, well-known material of single-crystal Bi will be helpful in the further research of topological materials.

FIG. 1 :
FIG. 1: (a,b) Angle-resolved photoelectron spectroscopy intensity plots along Γ-M taken at 7.5 K. Dashed lines indicate the edge of projected bulk bands.These figures are reproduced from Ref. 17. (c) Schematic drawing of the three-dimensional Brillouin zone (solid line) of the Bi single crystal and its projection onto the (111) surface Brillouin zone (dashed line).(d-g) Schematic drawings of the surface-state dispersion and projected bulk bands along Γ-M , obtained by these various methods: (d,e) Based on local density approximation with slab geometry including (c) and excluding (d) the interaction between the top and bottom surface.(f) Based on empirical tightbinding (TB) parameters and transfer-matrix method.(g) The same as (e) but using different TB parameters.
FIG. 2: Crystal structure of Bi.The lower contrast circles (light blue) represent the second atoms in the primitive unit cell (see text) in the (a) top and (b) side view of Bi(111).
(3.1)  to(3.3)  in Ref.16 , together with the appendix in Ref.15 .Any bulk states are the eigenstates of the 16×16 TM with unimodular eigenvalues.For each E in the projected bulk bandgap, T (k , E) has eight eigenvalues with moduli larger than 1, which correspond to the electronic states whose amplitude decays in the −z direction.With the boundary condition φ 0,1 = 0, the surface states should also decay outside the crystal.These surface states are determined by forming an 8×8 matrix M(k , E) composed of the eight components of φ 0,1 for each of the eight decaying states.The detailed procedure to generate M(k , E) is shown in Ref.22.Finally, the surface-state band dispersion (E(k )) is determined by solving det[M(k , E)] =0.

Figure 3 FIG. 3 :
Figure 3 plots the surface states and the edge of the projected bulk bands generated using the TM method with three different sets of bulk TB parameters.To plot the surface-state dispersion, we plotted ln(1/det[M(k , E)]), so that the locations (k , E) where surface states lie possess much smaller intensity values (i.e., darker) than the others .The surface bands calculated with the TB parameters given in Ref. 15 (Figs.3(a) and 3(b)) agree with those given in the previous paper exhibiting an artificial crossing of the surface bands between Γ and M owing to an incorrect mirror Chern number [16].This artificial surface-band crossing disappears when the new TB parameters are used.Using the parameters that generate topologically trivial bulk bands (TBP-1 in TableI), the surface-

FIG. 4 :
FIG. 4: (a) Band evolution at L. Solid (dashed) lines are the energies of bulk bands just above and below the bandgap at L calculated using the tight-binding (TB) parameter set TBP-1 (TBP-2) in Table I that give topologically non-trivial (trivial) topological order .Vertical lines indicate the position where the topological phase transition occurs, corresponding with each TB parameter set.(b-e) Surface-state band dispersions calculated by the transfer-matrix method with in-plane lattice distortions of (b,d) −5 % lattice strain and (c,e) +5 % lattice expansion.The dispersions are obtained using the TB parameter set (b,c) TBP-1 and (d,e) TBP-2, given in TableI.

FIG. 5 :
FIG. 5: Band dispersions calculated for a 20 biatomic-layer (BL) slab using the new tight-binding parameters that generate the (a) trivial (TBP-1 in TableI) and (b) non-trivial (TBP-2 in TableI) bulk topological order.Dashed lines represent the edge of the projected bulk bands.Intensities are obtained from the eigenfunction amplitude localized in the topmost surface BL.

FIG. 6 :
FIG. 6: The electronic structure around M in a 200 biatomic-layer Bi slab using the new tightbinding parameters that give (a) trivial (TBP-1 in Table I) and (b) non-trivial (TBP-2 in Table I) bulk topological order (cf.Figs.5(a) and 5(b), respectively).(c) Energy evolutions of the quantum-well-like states at M connected to the surface-state bands in the projected bulk band gap.Thick horizontal lines indicate the energies of the bottom of the projected bulk conduction band and the top of the projected bulk valence band.
Table I represents two sets of the new TB parameters, labeled as TBP-1 and TBP-2 in the following, obtained by the

TABLE I :
The new sets of tight-binding (TB) parameters for single-crystal Bi tuned so that the they could reproduce the surface states which agree with the experimental results.The definitions of each parameter are the same as in Ref.15.
methods above.TableIIrepresents the parity invariants at each TRIM in the Brillouin zone of the bulk Bi crystal, calculated with the parameters given in TableI.As shown the ν 0 values in Table