Analytical Solutions for the Surface States of Bi$_{1-x}$Sb$_x$ ($0\le x \lesssim 0.1$)

Analytical solutions for the surface state (SS) of an extended Wolff Hamiltonian, which is a common Hamiltonian for strongly spin-orbit coupled systems, are obtained both for semi-infinite and finite-thickness boundary conditions. For the semi-infinite system, there are three types of SS solutions: (I-a) linearly crossing SSs in the direct bulk band gap, (I-b) SSs with linear dispersions entering the bulk conduction or valence bands away from the band edge, and (II) SSs with nearly flat dispersions entering the bulk state at the band edge. For the finite-thickness system, a gap opens in the SS of solution I-a. Numerical solutions for the SS are also obtained based on the tight-binding model of Liu and Allen [Phys. Rev. B, 52, 1566 (1995)] for Bi$_{1-x}$Sb$_x$ ($0\le x \le 0.1$). A perfect correspondence between the analytic and numerical solutions is obtained around the $\bar{M}$ point including their thickness dependence. This is the first time that the character of the SS numerically obtained is identified with the help of analytical solutions. The size of the gap for I-a SS can be larger than that of bulk band gap even for a"thick"films ($\lesssim 200$ bilayers $\simeq 80$ nm) of pure bismuth. Consequently, in such a film of Bi$_{1-x}$Sb$_x$, there is no apparent change in the SSs through the band inversion at $x\simeq 0.04$, even though the nature of the SS is changed from solution I-a to I-b. Based on our theoretical results, the experimental results on the SS of Bi$_{1-x}$Sb$_x$ ($0\le x \lesssim 0.1$) are discussed.


Introduction
Spin-orbit coupling (SOC) produces an unconventional surface state (SS) that is inseparable from the electronic states of the bulk solid the surface belongs to. This unconventional SS is very different from SSs that are separable from the bulk states, which can be due to dangling bonds and surface reconstruction. From the viewpoint of the SOC-originated SS, bismuth is one of the most fascinating materials [1][2][3][4][5][6][7][8][9] because it has the largest SOC among the radioactively safe elements. Moreover, following the proposal of the existence of a threedimensional topological insulator in Bi 1−x Sb x , 10) the SS of Bi 1−x Sb x (especially for 0 ≤ x 0.1) has attracted considerable attention in solid-state physics. [11][12][13][14][15][16][17][18][19][20][21][22][23] While it has been widely accepted that pure Bi is topologically trivial, it has been proposed that Bi 1−x Sb x becomes topologically nontrivial 10,12) due to band inversion at the L point in a bulk Brillouin zone. [24][25][26][27][28][29][30][31][32] A drastic change in the SS is therefore expected due to this topological transition caused by the substitution of Sb into Bi.
The (111) surfaces of Bi 1−x Sb x have been examined extensively. While it seems that the experimental results for SS are heading toward to be settled, the interpretation of the results is still controversial; the SS around theM point is of a particular interest. (TheM point in a surface Brillouin zone corresponds to the L point in a bulk Brillouin zone.) According to the previously obtained experimental results, [1][2][3][4][5][6][7][8][9][11][12][13][14][15][16][17][18][19][20][21][22][23] there are two SSs between theΓ andM points in a surface Brillouin zone. Near theM point, one SS enters into a bulk conduction band, while the other enters into the bulk valence band. The controversy centers around the fact that the SS does not exhibit a qualitative change with respect to the Sb substitution for 0 ≤ x 0.2, even though the bulk conduction and valence bands at the L point are inverted at x ≃ 0.04, [24][25][26][27][28][29][30][31][32] where the topology of the system is changed. 10,12) (It should * fuseya@uec.ac.jp  be noted that in early studies on topological insulators, a third SS seemed to appear by the band inversion, 11,13,[16][17][18] however, recent experiments have attributed this third band to surface imperfections. 14,20) ) Because of these unexpected behaviors, it has recently been suggested that pure Bi should be topologically nontrivial. 19,22,23) If pure Bi is topologically nontrivial, then Bi 1−x Sb x should become topologically trivial after the band inversion (x 0.04), namely, a qualitative change in the SS should occur through the band inversion after all, even if pure Bi is topologically nontrivial. However, the experimental observa-tions that have been made do not indicate such a qualitative change. [15][16][17][18]20) Furthermore, the band structure of pure Bi was meticulously re-examined using first-principles calculations, and pure Bi was nevertheless determined to be topologically trivial. 21,33) The core of the mystery is at the properties of SS around theM point. Thus the problem should be solved by detailed investigation of the local region in the Brillouin zone. In this paper, we throughly investigate the SSs around theM point on the basis of k· p theory, an approach that has never before been adopted for the current problem. We succeeded in obtaining analytical solutions for the effective Hamiltonian common to the strongly spin-orbit coupled system -the Wolff Hamiltonian, the validity of which has been proved by the long history of studies on Bi and Bi-Sb alloys. [34][35][36][37] This approach is unbiased against the topological properties of the system. [38][39][40][41][42][43][44][45] The essence of our findings is as follows. For a semi-infinite system (i.e., the bulk limit), there are two solutions, I-a and I-b, which are responsible to the current the problem. In solution I-a, the two linear SS dispersions cross over one another in the direct bulk band gap; this is depicted in Fig. 1 (c). In solution I-b, on the other hand, one SS linearly enters into the bulk conduction band away from the band-edge, while the other enters into the bulk valence band as in Fig. 1 (f). This result has not been pointed out by the previous analytical approaches [38][39][40][41][42][43][44][45] and this study is the first to show it. Their existence condition is given by the parameter β = α ′ zz /α zz , where α zz is an inverse mass tensor that originates from the matrix elements between the bulk conduction and valence bands perpendicular to the surface, and α ′ zz is that from the other bands. Solution I-a is obtained when β < 0, while I-b is obtained when β > 0. For a system of finite-thickness, however, a gap opens in solution I-a, as can be seen in Fig. 1 (b); this gap can exceed the bulk band-gap as is shown in Fig. 1 (a); however, the SS of solution I-b does not change so much as the thickness is reduced as in Fig. 1 (d)-(f). Consequently, there is no apparent change in the SS even if the nature of SS is changed from I-a to I-b for systems with a finite-thickness [cf. Fig.  1 (a) and (b)]. These analytic results are seen to be in perfect agreements with the numerical simulation based on the model of Bi and Sb by Liu and Allen. 46) The remainder of this paper is organized as follows: in Section 2, the effective Hamiltonian for a strongly spin-orbit coupled system is introduced, and a differential equation with boundary conditions for the SS is derived; in Section 3, the results for a semi-infinite boundary condition are shown, while those for a system with finite-thickness are given in Secion 4 -up to this point, the analytical solution have been general so they can be applied to a variety of systems with strong SOCs; in Section 5, we apply our analytical results to Bi 1−x Sb x (0 ≤ x 0.1) and compare them with numerical results; Section 6 discusses the experimentally observed SSs of Bi 1−x Sb x based on our results; and Section 7 concludes this paper.

Model
It is known that the effective Hamiltonian for a bulk system with strong SOC is generally given by the Wolff Hamiltonian, 34,35) which is essentially to the same as the Dirac Hamiltonian but with spatial anisotropy. 37) The Wolff Hamiltonian is obtained by applying k · p theory only to the conduction  Illustration of the boundary conditions and the probability density |ψ(z)| 2 for (a) a semi-infinite system and (b) a system with a finite thickness. In this figure, the z-axis is considered to be perpendicular to the surface. and valence bands of a system, i.e., a two-band model. With this approximation, physical quantities such as the diamagnetic orbital susceptibility of Bi can be quantitatively calculated. 36,37) However, contributions from other bands, which are neglected in the original Wolff Hamiltonian, can in some instances play crucial roles. 47) The SS that originates from the SOC is one such example. In order to take into account the contributions from the other bands, in the present study we extend the Wolff Hamiltonian by using the Löwdin partitioning up to the second order. 48,49) The resulting "extended" Wolff Hamiltonian can be written as in the following form is the contribution from the other bands, while α ′ is an inverse mass tensor, and σ µ is the Pauli matrix. The contributions from the other bands are not generally symmetric between the conduction and valence bands. In this article, however, we begin with the above symmetric extended Wolff Hamiltonian, which is so general that it combines existing effective models for topological insulators and strongly spin-orbit coupled systems, such as HgTe quantum wells, 50) Bi 2 Se 3 , 42,51) and SnTe. [52][53][54] However, because of this generality, a more mathematically complex procedure is required. For example, in the previous models, −i ∇ z appears only in the σ z term, but it also appears in the σ x and σ y terms in our extended Wolff Hamiltonian. We analytically investigated the bound states near the surface for the Hamiltonian (1). A similar approach has been used for topological insulators, such as the HgTe/CdTe quan-tum wells 38,39) and Bi 2 Se 3 family [40][41][42][43] with each effective Hamiltonian. These works only considered solutions for the topologically nontrivial SS, i.e., a SS with a massless Diraclike dispersion. In the present study, the differential equation for the SS is solved without any assumptions being made about the topological nature of the system, i.e., both the topologically trivial and nontrivial SSs are treated equally, which allows for new SS solutions to be found.
The z-axis was set so as to be perpendicular to the surface, and we consider two situations: (i) a semi-infinite system in which the region z > 0 is the bulk solid under consideration, and z < 0 is a vacuum [ Fig. 2 (a)]; and (ii) a thin film with a finite thickness L, where |z| ≤ L/2 is the thin film, and |z| > L/2 is a vacuum [ Fig. 2 (b)]. The Hamiltonian for these systems is commonly given by the replacement p z → −i ∇ z . The only difference between these two situations is that of the boundary conditions used. The in-plane wave vectors k x and k y are good quantum numbers. Hereafter, p x,y will not be considered to be an operator, but rather an abbreviation of p x,y = k x,y . By introducing an ansatz wave function of the form ψ ∝ e λz into the extended Wolff Hamiltonian, Eq. (1) becomes the following form where λ = Λ/ (> 0) is the inverse of the localization length, and (The subscripts µ are assumed to be summed over pairs of identical indices.) By separating ǫ ′ ( p) into parallel (x, y) and perpendicular (z) components of the plane, we get the following equations In these equations, tilting in the z-direction is neglected. The eigenvalues of this Hamiltonian are easily obtained by considering its square, H 2 ψ = E 2 ψ. By using the relation where the inverse mass and the in-plane dispersion for the conduction and valence bands are given by The energy of the bulk state is given as |E| ≥ E bulk , where the following equation applies From Eq. (7), we can obtain the inverse localization length as a function of E in the form where the sign of the square root is defined to as + for Λ 1 and − for Λ 2 . Hereafter, the unit of the energy is set to be ∆ = 1.
Because the eigenvalues are doubly degenerate, the corresponding wave function can be written as a linear combination of two linearly independent eigenfunctions: where By considering the degrees of freedom of Λ 1,2 and ± signs in e ±Λ n z/ , the general solution for the wave function is given by a linear combination of eight eigenfunctions as 43) Ψ(z) = l=1,2 m=± n=1,2 C lmn ψ lmn e mλ n z .

Semi-infinite system
In this section, we focus on the solutions for the semiinfinite system. The boundary conditions for the semi-infinite system are given by the following equation: The latter condition requires m = −, so that the wave function is practically given as a linear combination of four eigenfunctions The condition Ψ(z = 0) = 0 yields the four simultaneous equations for C l−n , This leads to the eigenfunctions in the form The depth dependence of the probability density |ψ(z)| 2 of the wave function is illustrated in Fig. 2 (a). The wave function is localized very close to the surface and decays exponentially in the bulk region. The peak position depends on the difference between λ 1 and λ 2 . The secular equation of Eq. (21) is the following: where I is a 2 × 2 unit matrix. From this secular equation, we can obtain the following equation with respect to E: By substituting Eq. (11) into Eq. (24), we obtain an equation to solve in the following form which characterizes the relative curvature of the other band contributions to that of the conduction and valence bands in the direction of the z-axis. Its sign is closely related to the topological property of the model, 45) but there is no exact relationship between them. The solutions of Eq. (25) are as follows (The solution of E s2+ is obtained by using two linearly independent eigenfunctions similar to Eqs. (12) and (13) but

Solution I
There are two types of solution I. Their existence conditions depend on the sign of β as (I-a): Although both solution I-a and I-b have the same dispersion of ± 2ξ, which is linear with k, the regions in which they appear are different, as shown in Fig. 3. In the following figures, the in-plane dispersion is set to be ξ = α 2 k 2 /2 and ξ ′ = α ′ 2 k 2 /2, which does not lose any generalities. The wave number and the localization length are scaled in terms of ℓ 0 = / √ m∆ (in the case of Bi with 2∆ = 13.6meV, 46) ℓ 0 is estimated to be ℓ 0 = 33.3Å, which is about eight layers. The thickness of one Bi bilayer is 3.94Å. 1) ) Solution I-a appears around k = 0, and the linear dispersions cross at k = 0, i.e., they cross as the massless Dirac dispersion. Note that the origin of k is not necessarily at the center of the Brillouin zone; it can be any extremum of the band as the conventional k · p theory.
On the other hand, solution I-b can only appear away from  k = 0 since 1 + ξ ′ < 0 is never satisfied near k = 0. The existence of solution I-b has not been mentioned in previous works; [38][39][40][41][42][43][44][45] it has been recognized that no SS appears for β > 0. The existence of solution I-b is shown for the first time in the present work.
The important quantity of a SS is its localization length. We can calculate the inverse of the localization length Λ by substituting the solution into Eq. (11). The inverse localization length of solution I is obtained as The dispersion of solution I enters into the bulk state at ξ ′ = −1. At this k-point, Λ 2 becomes zero, and the SS naturally continues to the bulk state. It is interesting to note that the energy dispersion depends only on ξ, while the localization length depends only on ξ ′ . The probability of the wave function ψ(z) ≡ e −Λ 2 z/ −e −Λ 1 z/ of solution I-a is shown in Fig. 4 for different values of kℓ 0 in the case where β = −0.1, α = 1.0/m, α ′ = −0.5/m. The wave function is mostly bound at z/ℓ 0 ≃0.1-0.2, i.e., the first and second layers in the case of Bi. In this case, the SS enters into the bulk state at kℓ 0 = 2, where the localization length is infinitely long (Fig. 4), i.e., it naturally continues to the bulk  wave function. The properties of |ψ(z)| 2 for solution I-b are essentially the same as those for solution I-a.

Solution II
The energy dispersion of solution II is tangential to E bulk at the band edge E = ±1 and k x,y = 0. When ξ ′ − βξ = 0, the surface band is perfectly flat. When ξ ′ − βξ 0, the surface band has a finite mass that depends on the value of ξ ′ − βξ; it can be positive [ Fig. 5 (a)] or negative [ Fig. 5 (b)].
For solution II, the inverse localization length can be obtained as In the limit of k → 0 (ξ, ξ ′ → 0), where E s2 enters into the bulk state, Λ 2 becomes zero, and the localization length becomes infinitely long. The probability of the wave function, |ψ(z)| 2 , of solution II is shown in Fig. 6 as a function of the depth z. The properties of |ψ(z)| 2 are essentially common to the properties obtained by another set of parameters.
Note that the solution II with a nearly flat dispersion resembles the SS (or the edge state) with flat band in the graphene ribbons 55,56) or in the "line node semimetals" by Burkov- Hook-Balents. 57) They share a common feature that the SS enters into the bulk state at the band edge. On the other hand, there is no specific condition for the appearance of solution II, while there is a certain condition for that in line node semimetals.

Finite-thickness system
In this section, we consider a finite-thickness system. The boundary conditions for a system with thickness L is given by the following equation: The wave function is given as a linear combination of eight eigenfunctions (18). For example, the terms including C 1+1 and C 1−1 are explicitly written as where Similarly, the terms including C 2+1 and C 2−1 are given as The rest of the terms can be obtained by simply changing n = 1 → 2 and their coefficients D 1∼4 → D 5∼8 correspondingly in the above form. The boundary condition yields the eight simultaneous equations for D 1∼8 and its secular equation (8 × 8) is given by the following equations: where E n = E + 1 + ξ ′ − α ′ zz 2 Λ 2 n , C n = cosh(λ n L/2), and S n = sinh(λ n L/2). This leads to the z-dependence of the eigenfunction as where F 1 and F 2 are coefficients of four components. This functional form satisfies the boundary conditions given in Eq. (32) (an example of this functional form is illustrated in Fig.  2 (b)). The determinant of the matrix L is given by Here, we used the following relations: One possible difficulty in utilizing the Wolff Hamiltonian is that the parameters W(µ), P µ , and Z µ are complex in general, which may cause a phase dependent result in numerical calculations. By using the above relations, however, we were able to succeed in writing the equation only in terms of the mass tensors, which are not affected by the phase of the wave functions. This is practically important for the numerical calculations. Following this, we can obtain the equation with respect to E:  (45) where T n = tanh(λ n L/2). In the limit of L → +∞, where T 1 , T 2 → 1, Eq. (45) becomes equivalent to Eq. (24); the solutions of Eq. (45) should therefore naturally continues to those of the semi-infinite system. Figure 7 shows the dispersion of solution I-a and I-b for different thicknesses L/ℓ 0 = 12, 15, 20, and 100. We can see that the linearly crossing SSs of I-a open the gap as the thickness decreases; this is due to the interference between the opposite sides of the surfaces, which is also found in other systems. 39-41, 45, 58) Figure 8 shows the α zz dependence of the energy of I-a at k = 0, indicating the size of the gap, E g . Interestingly, it is found that the size of the gap is scaled with respect to √ α zz /L as shown in Fig. 8 (b), i.e., it is independent from α ′ zz . As α zz increases (the effective mass decreases) and L decreases, E g increases. The gap of the SS becomes as large as the size of the bulk band gap, E g = 2∆, when √ α zz m = L/ℓ 0 . Solution I-b also exhibits a thickness dependence, though its effect is not so remarkable. The energy of the point at which the SS enters the bulk band increases as the system becomes thinner, as shown in the inset of Fig. 7 (b).

Solution I
The localization length is almost the same as that obtained for the semi-infinite system, since the E-dependence of Λ n , Eq. (11), is common between the semi-infinite and finitethickness systems, though there is a slight modulation due to the opening of the gap.

Solution II
We mentioned that the equation with respect to E for the finite-thickness system, Eq. (45), naturally becomes that of the semi-infinite one in the limit of L → +∞. If we take the opposite limit, L → 0, Eq. (45) becomes the following equaion: To be more precise, solution II is also a solution in the zerothickness limit. This indicates that there is no quantum size effect for solution II, which is similar to the edge state of graphene ribbons with a flat band. 55) However, such a solution does not appear in the numerical calculation in the next section, and does not seem to appear in experiments. We will therefore only concentrate on solution I-a and I-b hereafter.

Comparison with the numerical simulations
In this section, we will compare our analytical solutions with the numerical solutions using the tight-binding model  Evolution of the conduction and valence bands at the L point with respect to the Sb concentration, x, based on the VCA. The theoretical value of x is scaled up by a factor of 2 so that the comparison with the results from experiments can be done in a more straightforward manner. 47) of Liu and Allen 46) for Bi 1−x Sb x (0 ≤ x ≤ 0.1). Figure 9 shows the obtained band structure and its symmetry at the L point. The order of the band, from its bottom to its top, is {L s , L a , L a , L s , L s , L a }. This order is classified to be topologically "trivial", according to Fu and Kane. 10) For the substitution of Sb, we adopt a simple virtual crystal approximation (VCA), in which the tight-binding parameters are obtained by a linear extrapolation between pure Bi and Sb. The band evolution obtained in this manner agrees well with the results of previous experimental investigations. [24][25][26][27][28][29][30][31][32] The band inversion at the L point occurs at x = 0.043 according to the present calculation -this can be seen in Fig. 9 (b) -and its topology changes from trivial to nontrivial at this point. 10) It should be noted that the Sb content here x is scaled up by a factor of 2 in order for comparisons with the experimental results to be more straightforward. The details of the calculations are described in a previous study. 47) Figure 10 shows the energy dispersion of the finite bi-  100, 300, and 500 bilayer (BL) systems, as calculated using the Liu-Allen model. 46,47) The shaded region corresponds to the bulk band. The origin of the energy is fixed to being that of the Fermi energy of pure Bi for both (a) and (b). Note that pure Bi is topologically trivial and Bi 0.9 Sb 0.1 is topologically nontrivial for the Liu-Allen model according to the classification of Fu and Kane. 10) layer (BL) system around theM point based on the Liu-Allen model. For both Bi and Bi 0.9 Sb 0.1 , there are two SSs around theM point; however, their thickness dependences are different. For pure Bi, one of the SSs enters into the bulk conduction band, while the other enters into the valence band of a 100 BL system. In 300 and 500 BL systems, however, the SSs do not enter the bulk band and tend to cross with one another in the direct gap of the bulk band. This is a characteristic feature of solution I-a. The thickness dependence of the energy at k = 0 for a finite BL system agrees perfectly with the analytical solution we obtained using Eq. (45) at k = 0 (which is shown in Fig. 11). For the analyticalal solution, the input parameters are only α zz = 498.7/m and α ′ zz = −3.554/m, which are obtained from the "bulk" calculation of the Liu-Allen model. The perfect agreement between them indicates the correctness and validity of the analytical solution. This is the first time that the nature of the SS in the numerical calculation is clarified with the help of analytical solusions. It is rather surprising that even the 500 BL system exhibits a gap due to the interference between both side surfaces, i.e., 500 BL (≃ 200 nm thickness) cannot be regarded as the bulk system. 23) This remarkable quantum size effect is due to the small bulk bandgap and the small effective mass of Bi.
For Bi 0.9 Sb 0.1 , however, the SSs enter into the bulk bands Liu-Allen model Fig. 12. Sb content dependence of β = α ′ zz /α zz obtained by using the Liu-Allen model with a VCA. even in the 500 BL system. They are considered to be solution I-b; if so, the sigh of β of the "bulk" Bi 1−x Sb x should change from negative to positive due to the band inversion. The Sb content dependence of β, which is based on the Liu-Allen model, is plotted as a solid line in Fig. 12. β can be seen in this figure to change its sign from negative to positive, which corresponds to the qualitative change in SS from solution I-a to I-b. This is a clear evidence for the qualitative transition of the SS due to the band inversion. Therefore, as far as based on the Liu-Allen model, the sign change in β accompanies the band inversion.
The largest contribution to β for the conduction band at the L point in pure Bi is that of the L a band located at 0.9 eV shown in Fig. 9 (a). It may be rather surprising that the contributions from such a high energy band more than 1 eV far from the Fermi energy drastically affect the appearance SS. However, this can happen in the strongly spin-orbit coupled system, in which the large SOC (∼ 1.8 eV in Bi) produces a significant interband effect. 47) The most important consequence is that the appearance of the SS in the finite-thickness system is essentially unchanged by the band inversion even if the true nature of the SS is changed. For example, when we compare the 100 BL system results of pure Bi and Bi 0.9 Sb 0.1 in Fig. 10, we cannot dis-tinguish between them by only looking at their dispersions; rather, the differences only appear in very thick BL or semiinfinite systems. According to the thickness dependence exhibited in Fig. 11, the difference can be seen in the system with more than 200 BL (≃ 80 nm) at least.

Comparison with the experimental results
First, we briefly summarize the experimental results on Bi 1−x Sb x (0 ≤ x 0.1). The bulk band inversion in Bi 1−x Sb x has been examined repeatedly more than five decades through various types of measurements, such as magnetooptics, 26,30,31) transport and quantum oscilations, 24,25,27,29) and magnetic susceptibility. 28 The finite-thickness 8,9,18,20,22) and the semi-infinite [4][5][6]19) systems exhibit these properties in common. The essence of these experiments is the fact that the SSs seem to remain qualitatively unchanged through the band inversion, at which the topology of the system is changed. The results of present analytical studies supported by numerical calculations are in complete agreement with these experimental indications.
As has been discussed in the previous sections, when the thickness is finite, a gap opens in solution I-a. In the case of Bi 1−x Sb x , because the bulk band-gap and the bulk band-mass are very small, this quantum size effect appears to be drastic enough that the gap of the SS exceeds the gap of the bulkband. In reality, the gap of the SS in a 200 BL system is comparable to the bulk band gap, according to the estimation we made. (Note that the most recent measurements were carried out up to 202 BL for pure Bi 22) and up to about 300 BL for x > 0.07. 20) ) Therefore, for the system with finite thickness ( 200), solution I-a for x 0.04 cannot be distinguished from solution I-b for x 0.04. Therefore, the properties of SSs are apparently unchanged by the band inversion.
The difference between solution I-a and I-b should appear in the semi-infinite system. When the band inversion is accompanied by a sign change in β, there should be a qualitative transition from solution I-a to I-b in case of the bulk surface. If β < 0 in pure Bi, as is estimated by the Liu-Allen model, then linearly crossing SSs should be observed before the band inversion (0 ≤ x 0.04). However, such a gapless SS has not yet been observed, so there still remains an inconsistency between our numerical simulations based on the Liu-Allen model and the experiments on the bulk surface.
This discrepancy between our theoretical results and the experimental results for the bulk SS can be resolved if the resolution of the surface ARPES is improved. In pure Bi, the bulk band gap is 11-15 meV, and it is much smaller in Bi 1−x Sb x .
On the other hand, the current experimental resolution of surface ARPES was 25 meV in the study conducted by Ast and Höchst, 5,6) and it was 7 meV in the study conducted by Ohtsubo et al., 19) which are not high enough for Bi 1−x Sb x .
It is worth noting the recent progress made for ARPES measurements of the SS of the iron-based superconductor FeTe 1−x Se x . 59) Although the first-principles calculations predict a bulk band-gap of ∼ 10 meV and the linearly crossing SSs, such a SS has not been identified by ARPES measurements that had energy resolutions of ∼ 5 meV and ∼ 10 meV. 60) However, by using high energy and momentum resolution ARPES with an energy resolution of ∼ 70 µeV, linearly crossing SSs and the bulk band-gap were clearly detected. 59) Similar progress is therefore expected for Bi 1−x Sb x .

Some remarks on the topology of the system
To determine the topology of the system, we need to take into account the entire Brillouin zone; especially the parity eigenvalue at eight time-reversal invariant momenta. On the other hand, since our theory is based on k · p theory, which is rigorous for the local region in the Brillouin zone, we cannot give a statement on the topology of the system. Nevertheless, we can make some remarks on the topology of the system based on our results as follows.
In some cases, the topology of the system has been determined only by counting the intersections of the SSs and the Fermi energy. However, argument of such kind can lead to an incorrect solution for the system with finite thickness, since the opening the gap of SS due to the interference between opposite sides of surfaces can change the number of intersections, even though the topology of the system is unchanged. In connection with this, the conclusion that the pure Bi must be topologically nontrivial since one of the two SSs enters into the bulk conduction band and the other enters into the bulk valence band near theM can be incorrect for the system with finite thickness. There is a possibility that even topologically trivial system exhibits such SSs entering conduction and valence bands separately as is shown in Fig. 10 (a).

Conclusion
We have succeeded in obtaining the exact solutions for SSs of both semi-infinite and finite-thickness systems described by the extended Wolff Hamiltonian, which is common to strongly spin-orbit coupled systems. There are three types of solutions for the semi-infinite system; solution I-a and Ib have a dispersion linear in k, and solution II has a nearly flat dispersion.
In the semi-infinite system, solution I-a appears when β < 0 and 1 + ξ ′ > 0. The dispersion of solution I-a linearly cross with each other at k = 0, where k is measured from an extremum of the bulk band. On the other hand, solution I-b appears when β > 0 and 1+ξ ′ < 0. One of the SSs of I-b linearly enters into the bulk conduction band, and the other enters into the bulk valence band away from k = 0. Solution I-a corresponds to the gapless Dirac-like SS that has previously been discussed for topological insulators. Solution I-b is a finding that has been made for the first time. The key parameter that classifies the nature of a SS (solution I-a or I-b) is the sign of β.
In a finite-thickness system, a gap opens in solution I-a due to the interference between opposite sides of the surface. The size of the gap is scaled by √ α zz /L, and it becomes comparable to the bulk band gap for √ α zz m = L/ℓ 0 . As a result, this quantum size effect is more significant for smaller effective masses and bulk band gaps. We have also calculated the SS of Bi 1−x Sb x for 0 ≤ x ≤ 0.1 based on the tight-biding model of Liu and Allen. We have obtained a perfect agreement between the analytical and numerical results: the thickness dependence of the gap of solution I-a, and the qualitative change in the SS from solution I-a to I-b through the band inversion at x ≃ 0.04. Due to this agreement, the nature of the numerically obtained SS has been identified for the first time based on the analytic solutions. For a system with a thickness below about 200 BL (∼80 nm), the gap in solution I-a exceeds the bulk band gap. As a result, solution I-a becomes qualitatively the same as solution I-b. This can give a possible answer to the experimental mystery on the (111) surface of Bi 1−x Sb x for the finite-thickness system. However, there still remains discrepancy between our results and experiments for the semi-infinite system. For this discrepancy, high-resolution ARPES measurements will give a crucial progress.
Although we focused on the problem on Bi 1−x Sb x in this paper, the present analytical solutions are not restricted to this specific case. They are applicable to various systems with strong SOC, including IV-VI semiconductors (PbTe, SbSe, SnTe), and topological insulators.