Momentum-space structure of surface states in a topological semimetal with a nexus point of Dirac lines

Three-dimensional topological semimetals come in different variants, either containing Weyl points or Dirac lines. Here we describe a more complicated momentum-space topological defect where several separate Dirac lines connect with each other, forming a momentum-space equivalent of the real-space nexus considered before for Helium-3. Close to the nexus the Dirac lines exhibit a transition from type I to type II lines. We consider a general model of stacked honeycomb lattices with the symmetry of Bernal (AB) stacked graphite and show that the structural mirror symmetries in such systems protect the presence of the Dirac lines, and also naturally lead to the formation of the nexus. By the bulk-boundary correspondence of topological media, the presence of Dirac lines lead to the formation of drumhead surface states at the side surfaces of the system. We calculate the surface state spectrum, and especially illustrate the effect of the nexus on these states.

The possible topological defects occurring in gapless materials are not limited to these options. The next step in the increasing complexity of the momentum space topological defects is to consider the topology of the systems containing multiple Dirac lines 67,68 . These Dirac lines can meet somewhere in the momentum space giving rise to an exceptional point, which can be called as nexus 67 . In real space nexuses can appear when vortex lines meet and such kinds of topological defects have been considered in the contexts of superfluid Helium-3 1 and as a possible mechanism for confinement in the Georgi-Glashow model in particle physics 69 . Recently, the idea (b) Momentum-dependent energy gap between the two bands closest to the Fermi level E3(k)−E2(k) as a function of kx and ky for kz = 0. There are band crossings in the vicinity of K and K points in the (kx, ky)-plane (Fig. 3). As a function of kz these band crossings form nodal Dirac lines, so that for each momentum along the line the energy spectrum with respect to the transverse momenta is conical. The band crossings always appear within the mirror planes directed along the kz direction and the dashed lines within the (kx, ky)-plane, and are protected by structural mirror symmetries. The mirror planes share common lines in the momentum space directed along kz-direction exactly at K and K in the (kx, ky)-plane. The Dirac lines can merge at a specific kz within this line, leading to an appearance of a nexus.
that a momentum space nexus may occur in topological semimetals was put forward in Ref. 67 in the context of Bernally (AB) stacked graphite, where the existence and merging of multiple band contact lines is theoretically well-established [70][71][72] .
In this paper we study the properties of the nexus semimetal phase proposed in Refs. 67, 70-72. We identify the key properties of the model that allow stabilizing the nexus in the momentum space: reflection symmetries and the existence of more than one electron-or hole-like bands close to the Fermi energy. The reflection symmetries of the structure specify mirror planes in the momentum space where the Hamiltonian is block-diagonal, and allow stabilizing band crossings where the bands have different eigenvalues of the mirror symmetry operator 16,48 . The requirement for stabilizing the nexus is to have several of those mirror planes sharing a common line in the momentum space, so that stable Dirac lines can merge at a specific point within this line. We show that in the nexus semimetal phase proposed in Refs. 67, 70-72 the nexus indeed defines a point-like topological defect in the momentum space where three bands are degenerate at the same momentum, and the appearance of this type of an exceptional point distinguishes the nexus semimetal phase from a semimetal containing (multiple) Dirac lines that do not merge in the momentum space. By the bulk-boundary correspondence of topological media, the presence of Dirac lines lead to the formation of surface states. We show that the number of surface states as a function of the momentum components parallel to the surface displays a fine-structure associated with the projected Dirac lines. Far away from the nexus the surface state dispersions take a form of a drumhead that is bounded by the projected crossing points of the electronand hole-like bands 16,48 . In the vicinity of the nexus there occurs a crossover to a different type of behavior, where one of the surface bands connects two electron-like bands to each other, and the other surface band connects two hole-like bands. In the crossover regime the surface states hybridize with the bulk states so that they connect bulk band edges to each other instead of being bounded by the projected Dirac lines.

II. MODEL AND SYMMETRIES
We consider Bernally stacked honeycomb lattices [ Fig. 1(a)] such as Bernal graphite 67,[70][71][72] . The tightbinding model for such kind of three-dimensional system [in the layer ⊗ sublattice space described in Fig. 1(a)] can be written as The different hopping parameters γ i are illustrated in Fig. 1(a). We have neglected further neigbor hoppings, as they do not change the qualitative conclusions if they are small and preserve the structural symmetry. Additionally, there exists a parameter ∆, denoting a locally broken A-B sublattice symmetry but still preserving the global A-B symmetry. This term is allowed by the symmetries of the structure, since one of the sublattices in each layer has an atom on top of it so that it has a different environment than the other sublattice. The structure factors arising from the Fourier transform are f (k x , k y ) = e i δi·(kx,ky) and Γ(k z ) = cos(k z ). The nearest neighbor vectors δ i [ Fig. 1(a)] are normalized so that the vectors connecting neighboring unit cells have unit length. In z-direction we use the spacing between the layers as the unit length so that −π/2 ≤ k z ≤ π/2.
The most important symmetries of the model (see App. A for more details) are a mirror symmetry and a three-fold rotational symmetry where (k x ,k y ) and (k x ,k y ) are the momentum coordinates after a rotation by ±2π/3 around the z-axis. There are similar mirror symmetries also with respect to (k x ,k y ) and (k x ,k y ). In a special case ∆ = γ 4 = 0 the system supports an accidental chiral symmetry CH(k)C = −H(k), where C = τ 0 σ z . In graphite, ∆, γ 4 γ 0 , γ 1 , γ 3 , so that the chiral symmetry is valid as a good approximation.
There are special planes going through the Γ point and at the boundary of the Brillouin zone which are mapped back to themselves in the mirror symmetries (up to a reciprocal lattice vector) [ Fig. 1(b)]. The relevant three planes around the K point are directed along the k zdirection and k y = 2π/ √ 3,k y = −2π/ √ 3 andk y = 0 within the (k x , k y )-plane 73 . Within these mirror planes the mirror symmetries give rise to symmetries commuting with the Hamiltonian at fixed momentum, e.g.
The symmetry operators in different coordinates are , and all of them are simultaneously valid within the line along k zdirection at the K point: In addition to the symmetries we assume a hierarchy of couplings |γ 0 | |γ 1 | > |γ 3 |, |γ 4 |, |∆|. In particular, this fixes the overall behavior of the different bands so that two of the bands are electron-like (hole-like), bending upwards (downwards) in energy when moving away from the K point. Some of the details discussed below depend on the relative signs of the couplings. In the main text we consider all couplings to be positive [in the convention defined by the Hamiltonian (1)] and we illustrate the main effects of other choices of signs in App. E. When not otherwise stated, we choose in the figures γ 1 = 0.3γ 0 , γ 3 = ∆ = 0.1γ 0 and γ 4 = 0.05γ 0 . These are close to the values 74 often used for Bernal graphite, except that we use somewhat larger ∆ and γ 4 to better illustrate the properties of the nexus.

III. BULK PROPERTIES
In the following we study the bulk spectrum around the K-point [(k x ,k y ) =K + (q x , q y )] 73 . BecauseSS = diag(1, e i2π/3 , e −i2π/3 , 1) andS commute with the Hamiltonian, two bands must always be degenerate at the K point for all values of k z . This degeneracy also follows from the tight-binding Hamiltonian (1) but it is protected by the existence of the mirror symmetries, not by the explicit form of the Hamiltonian. However, there are two topologically distinct options on how two bands can be degenerate within this model as illustrated in Fig. 2. For small k z electron and hole bands are degenerate at the K point [ Fig. 2(a)], whereas for large k z two hole bands are degenerate at the K point and there is an energy gap between the electron and hole bands [ Fig. 2(c)]. Because the mirror symmetries guarantee that at least two bands are always degenerate at the K point, there must exist a value of k z [Γ(k z ) = ∆/(2γ 1 )] where three bands are simultaneously degenerate [ Fig. 2(b)]. We call this exceptional point nexus because it is also a merging point of several Dirac lines (see below). The existence of the exceptional point where three bands are simultaneously degenerate is protected by the fact that situations (a) and (c) are topologically distinct in the presence of the mirror symmetries.
In addition to the band contact line at the K point, three other band crossings appear 67,[70][71][72] at energy E c within the three mirror planes at q ≡ (q x , q y ) = q c (−1, 0) and The two bands crossing at q = q c (−1, 0) correspond to different mirror eigenvaluesS|ψ ± = ±|ψ ± . There- Around the crossings at q = 0 and q = q c (−1, 0) the projected Hamiltonians in the basis of eigenvectors corresponding to eigenvalues +1 and −1 of the mirror symmetry operator at the crossings are 76 The latter describes a tilted Dirac cone for each value of k z , and it becomes overtilted at Γ(k * z ) ≈ ∆ γ1 . For k z > k * z , the Dirac line can thus be called type II Dirac line 77 (see App. A). These Hamiltonians are locally of the form where the direction of the pseudomagnetic field h 0[c] (δq) rotates by 2π when δq goes around the band crossing point [ Fig. 3(d)]. Therefore, the wave vector would obtain a Berry phase ±π if taken around such a path. Such Berry phases are typically associated with the appearance of surface states 48 . However, the Berry phase is only defined modulo 2π, and as a result this argument can be used to explain the existence of only a single protected surface state at a given surface momentum.

IV. SURFACE STATE SPECTRUM
We consider a translationally invariant system in the x-and z-directions, so that k x and k z are good quantum numbers. In the presence of chiral symmetry (∆ = γ 4 = 0) the 1D Hamiltonian H kx,kz (k y ) has a well-defined topological invariant (winding number) W (k x , k z ) whenever there is no gap closing as a function of k y (see App. B), and this invariant determines the number of zero-energy surface states for each k x and k z . Moreover, W (k x , k z ) can change only at the Dirac lines with where the energy gap closes 73 . By computing W (k x , k z ) in the presence of chiral symmetry (see App. B) we arrive at a flat band (zero energy) spectrum in the regions of the transverse momenta with W = 0 in Fig. 4(a). There are two flat bands connecting the projected Dirac lines In the presence of the chiral symmetry the Dirac lines meet at the boundary of the Brillouin zone k z = π/2 and there is no three-fold degenerate exceptional point where a gap between the electron and hole bands opens up. Therefore, ∆ and γ 4 have important effects on the surface state spectrum. The eigenstates in the presence of chiral symmetry can be solved exactly, and the surface state dispersions can be computed perturbatively in the limit γ 0 , γ 1 Γ(k z ) γ 3 Γ(k z ), γ 4 Γ(k z ), ∆. By utilizing also the observation that the surface bands connect the Dirac lines in energy and momentum we get (see App. C) for for K x − q c (k z ) < k x < K x , which are valid far away from the nexus 78 (Fig. 4). Therefore for small k z below the nexus the dispersions take the form of a drumhead that is bounded by the projected crossing points of the electron-and hole-like bands 16,48 [ Fig. 4(b)]. On the other hand, for large k z on the other side of the nexus there is a gap between electron-and hole-like bands and one of the surface bands connects two electron-like bands to each other, and the other surface band connects two hole-like bands. For k z = π/2 the surface state dispersions can be solved exactly for Hamiltonian (1) and one obtains E 1 = 0 and E 2 = ∆ for 2π/3 ≤ k x ≤ 4π/3 (see App. C) [ Fig. 4(d)]. This qualitative change in the behavior of the surface bands on opposite sides of the nexus signals the existence of an exceptional point in the momentum space where three bands are simultaneously degenerate. By numerically computing the surface Green functions 79 , we find that in the vicinity of the nexus the surface states hybridize with the bulk states so that they connect bulk band edges to each other instead of being bounded by the projected Dirac lines [ Fig. 4(c)]. This hybridization appears in the regime where the Dirac cones around the band crossings [Eq. (2)] are overtilted (|dξ c /dq x | > |d h c /dq x |) forming type II Dirac lines (see App. A and App. D) 77 .

V. SUMMARY AND DISCUSSION
We have identified the symmetries of the model that allow stabilizing the nexus in the momentum space and shown how the momentum-space structure of surface states follows from the properties of Dirac lines. In the vicinity of the nexus the behavior of the surface states changes qualitatively, indicating the existence of triple degeneracy point in the momentum space. There is an ongoing search for new types of fermions in condensed matter systems (in addition to Majorana, Weyl and Dirac fermions) that are described by simultaneous crossings of multiple bands 80 . The best candidate material for the study of nexus fermions is regular graphite, where the surfaces with a component parallel to the c-axis should exhibit surface states. The properties of these surface states can be studied with STM and ARPES. Other candidate materials include for example suitably stacked silicene layers 81 and InAs 1−x Sb x 82 . The latter also supports a pair of triple degeneracy points connected by Dirac lines. Moreover, this material obeys similar mirror and three-fold rotational symmetries as graphite, so that the stabilization of the nexus in InAs 1−x Sb x follows from our analysis. We point out that the symmetry analysis described in this paper may be a useful starting point for a general classification scheme of the nexus semimetals based on the space group symmetries, and it would be interesting to find out whether the predicted properties of the surface state spectrum in the vicinity of the nexus are generic for all nexus semimetals.
Apart from spectroscopic features, topological phases often have unusual response characteristics (anomalies) [83][84][85][86][87][88] . One interesting direction for future research is to find out whether the nexus semimetal phase is characterized by an anomaly associated with a spectral flow between the nexus points. We also point out that even relatively weak interactions can lead to symmetry-broken states at the surface because of the large density of states caused by the approximately flat bands [89][90][91][92][93][94] . These symmetry-broken states are expected to be exotic since they cannot be described with a mean field theory 95 . Finally, we expect that the interactions may lead to "dipole" correlations within the structure, so that the effective ∆ is renormalized. This would mean that the position of the nexus and the size of the energy gap above it would depend on temperature.
We consider the tight-binding model for Bernally Additionally, there exists (iv) a three-fold rotational symmetry , so that similar mirror symmetries exist also with respect to (k x ,k y ) and (k x ,k y ). In a special limit ∆ = γ 4 = 0 the system also supports a chiral symmetry CH(k)C = −H(k), where C = τ 0 σ z .
The most important symmetries are the mirror symmetries with nontrivial matrix structure H(k x , k y , k z ) = τ x σ x H(k x , −k y , k z )τ x σ x [and correspondingly for (k x ,k y ) and (k x ,k y )]. There are special planes going through the middle of the Brillouin zone and at the boundary of the Brillouin zone which are mapped back to themselves in the mirror symmetries (up to a reciprocal lattice vector). The relevant three planes around the K point are directed along the k z -direction and k y = 2π/ √ 3, k y = −2π/ √ 3 andk y = 0 within the (k x , k y )-plane. Within these mirror planes the mirror symmetries give rise to symmetries commuting with the Hamiltonian at fixed momentum The symmetry operators in different coordinates are By analyzing this Hamiltonian it is easy to see that there always exists a band crossing at (q x , q y ) = (0, 0) for all values of k z as discussed in the main text. Additionally there exists three other Dirac lines in the vicinity of K point. These band crossings appear at finite energy E c within the three distinct mirror planes at (q x , q y ) = q c (−1, 0) and (q x , q y ) = q c (1/2, ± √ 3/2), where and Around the crossings at q = 0 and q = q c (−1, 0) the projected Hamiltonians in the basis of eigenvectors corresponding to eigenvalues +1 and −1 of the mirror symmetry operator at the crossings are given by Eqs. (2) in the main text. The Hamiltonian around q = q c (−1, 0) describes a tilted anisotropic Dirac cone for each value of k z . For small values of k z the Dirac cone is only slightly tilted but when one approaches the nexus [Γ(k z ) = ∆/(2γ 1 )] the tilt increases. The Dirac cone becomes overtilted at For k z > k * z , the Dirac line can thus be called type II Dirac line (see also the discussion in Ref. 77). This transition, which occurs already at slightly smaller value of k z than the nexus, is illustrated in Fig. 5. We consider a translationally invariant system in the xand z-directions corresponding to a zigzag edge for each graphene layer. Therefore k x and k z are good quantum numbers and by fixing them we get a 1D Hamiltonian H kx,kz (k y ), which depends only on k y . In the special limit ∆ = γ 4 = 0, the system supports an additional chiral symmetry C = τ 0 σ z , and the 1D Hamiltonians H kx,kz (k y ) have well-defined topological invariants. To calculate this topological invariant, we first notice that the Hamiltonian can be written in a block-off-diagonal form where , (B2) = e 2i δ·(kx,ky) f (k x , k y ) and δ = (0, 1/(2 √ 3)). The topological invariant can then be defined as a winding number where the integration is over the 1D Brillouin zone in k y direction. The winding number undergoes a series of transitions at the momenta of the projected Dirac lines k These changes occur because the complex field z has vortex lines at the positions of the Dirac lines [ Fig. 6(a),(b)]. The winding number W as a function of k x and k z is shown in Fig. 4(a) in the main text.
The winding number W (k x , k z ) determines the number of zero-energy states for each k x and k z . In the presence of chiral symmetry, the band crossings always occur at energy E = 0, and the surface states form flat bands at E = 0 between the Dirac lines [ Fig. 6(c)].
In the absence of chiral symmetry the band crossings are protected by the mirror symmetry and they appear at finite energy E c . Moreover around some of the crossings the low energy theory is described by a tilted Dirac cone. By considering a general tilted Dirac cone we can introduce an edge by replacing q y = −i∂ y . By looking for an exponentially localized solution ψ = a 1 a 2 e −αy at the energy E s (q x ), we get Therefore real solutions of α exist only if and by varying E s (q x ) (and the sign of q x ) one can interpolate between different boundary conditions determining the ratio a 2 /a 1 . [The surface state may occur either for q x < 0 or q x > 0 depending on the boundary conditions. Similarly the exact dispersion E s (q x ) depends on the boundary conditions.] As one can see from inequalities (B7), for q x = 0 the surface state energy must satisfy E s (0) = E c , so that a single Dirac line gives rise to a drumhead surface state dispersion bounded by the projected Dirac line in energy and momentum. In particular, it follows from this calculation that the existence of a single surface state is always guaranteed in the regime Fig. 6(d)]. On the other hand, in the regime K x + q c (k z )/2 < |k x | < π we expect to find two surface states in the presence of chiral symmetry, and once the chiral symmetry is broken due to ∆, γ 4 = 0 the existence of the drumhead surface states may depend on how strongly the Dirac cones are tilted (see below).

Appendix C: Analytical results for the surface state spectrum
In order to obtain analytical insights into the surface state dispersions, we start by considering the zero energy wave functions in the presence of the chiral symmetry (γ 4 = ∆ = 0). The solutions exist either only in sublattice A or sublattice B depending on which surface one is considering. In the following we concentrate on those solutions which are localized in sublattice B where K y = 2π/ √ 3 is the y-component of the momentum at the Dirac point and the indices refer to the layer degree of freedom. By substituting the ansatz [Eq. (C1)] to the block off-diagonal form of the Hamiltonian [Eq. (B1)], we arrive at equations We look for a solution of the form (the plane y = y 0 describes the surface) The wave functions localized in sublattice B exist on the right surface, i.e., Re This solution describes how the weight of the wave function within the different layers varies as a function of k x and k z in the case of the flat band corresponding to W = 1 regime in Eq. (B4). Additionally we need to find expressions for the wave functions of the flat bands correponding to W = 2 regime. To describe these wave functions we look for solutions of Eq. (C3) allowing also complex values of Q, but still requiring Re[Q] < 0. Such solutions exist for K x + q c /2 < k x < π and they come in pairs Q and Q * , where The orthonormal solutions obtained using these solutions can be written as where In the special limit Γ(k z ) = 0, the Hamiltonian becomes block-diagonal in the layer degree of freedom, so that Q = Q * and these expressions describe two copies of the edge states for a single layer graphene. On the other hand for Γ(k z ) = 1 the solutions describe the edge states of bilayer graphene. By considering the couplings γ 4 and ∆ as a perturbation, we obtain in lowest order in γ 3 for K x + q c (k z )/2 < |k x | < π and Additionally, we can utilize the knowledge that far away from the nexus the surface states must have the shape of the drumhead bounded by the projected Dirac lines. By including the necessarily corrections to E 1 and E 2 , we arrive to surface states dispersions described by Eqs. (3) in the main text. In principle these expressions should work for Γ(k z ) ∆/2γ 1 or |k x − K x | ∆/γ 0 . However, numerically we find that they work everywhere except very close to the nexus. On the other hand, since the single protected surface state should exist also in the vicinity of the nexus, we find that the surface dispersion E 3 is given by Eq. (4) in the main text.
Finally we notice that the full tight-binding Hamiltonian also has an exact solution for Γ(k z ) = 0 i.e. for k z = π/2. Namely, in this case the Hamiltonian becomes block-diagonal in the layer degree of freedom, and each block can be solved similarly as the surface states for a single layer graphene. This way, we obtain E 1 = 0, E 2 = ∆ for Γ(k z ) = 0 and K x ≤ |k x | ≤ π.
(C10) These also coincide with Eqs. (3) for Γ(k z ) = 0. We have checked that the analytical solutions (3), (4) and (C10) describe the surface state dispersions reasonably far away from the nexus by numerically diagonalizing the tight-binding Hamiltonian in the case of a finite width in y-direction. However, in the vicinity of the nexus it is difficult to obtain analytic expressions for the surface state dispersions. Moreover, in that regime the localization length of the surface states becomes very long, and thus the numerical diagonalization of the tight-binding Hamiltonian with large enough width in y-direction also becomes computationally expensive.
Alternatively the surface state dispersions can be obtained by numerically calculating the surface Green func- (the matrix indices 00 correspond to the surface in the y-direction). The surface density of states is given by This method is computationally much more efficient since the Green function can be computed for a semi-infinite system using a quickly converging renormalization group method 79 . In the numerics we broaden the δ-peaks in the ρ s (E, k x , k z ) corresponding to the surface state energies E s (k x , k z ) by using a nonzero value of η. Sufficiently far away from the nexus, the Green function method reproduces the analytic surface state dispersions (3), (4) and (C10) as shown in Fig. 4(b) and (d) in the main text. For small k z the dispersions take a form of a drumhead that is bounded by the projected crossing points of the electron-and hole-like bands [ Fig. 4(b) in the main text], but for k z = π/2 one of the surface bands connects two electron-like bands to each other, and the other surface band connects two hole-like bands [Fig. 4(d) in the main text]. We now turn to the description of the transition between these qualitatively distinct regimes, which occurs in the vicinity of the nexus. Figure 7 shows the surface density of states ρ s (E, k x , k z ) as a function of k x and E for specific values of k z . The surface state energies E s (k x , k z ) show up as peaks in ρ s (E, k x , k z ) (broadened by η) and form lines as a function of k x . Sufficiently far from the nexus the surface state dispersions take a form of a drumhead that is bounded by the projected crossing points of the electron-and hole-like bands [ Fig. 7(a)]. When the nexus is approached by decreasing Γ(k z ) [Γ(k z ) → ∆/(2γ 1 )] the bulk Dirac cones become strongly overtilted and the surface states hybridize with bulk states so that they connect bulk band edges to each other instead of being bounded by the projected Dirac lines [ Fig. 7  for large k z two hole bands are degenerate at the K point and there is an energy gap between the electron and hole bands. In the case ∆ < 0, for small k z electron and hole bands are still degenerate at the K point [ Fig. 8(a)]. However, now two electron bands are degenerate at the K point for large k z [ Fig. 8(c)]. In general all the results remain qualitatively similar when the sign of ∆ is changed except that one needs to mirror all the figures in energy E → −E.
The sign of γ 4 mainly influences the band crossing energies E c (k z ) as can be seen from Eq. (A2). Since the drumhead surface states connect band crossings in energy and momentum, γ 4 also influences the surface state dispersions [Eqs. (3) and (4)]. In particular the relative sign of γ 4 and ∆ determines whether the drumhead dispersion in the regime K x − q c (k z ) < |k x | < K x is tilted in the same direction for all values of k z (unidirectional surface states). This can for example influence the transport properties of these systems.
The sign of γ 3 determines the side in which the Dirac lines are with respect to the K-point [cf. Fig. 3(a) in the main text and Fig. 9(a)]. As a result the winding number W (k x , k z ) in the presence of chiral symmetry (∆ = γ 4 = 0) is modified in the different regions of k x and k z around k x ≈ K x [cf. Fig. 4(a) in the main text and Fig. 9(b)]. As discussed in App. B, the winding number influences the number of surface states also when the chiral symmetry is broken (∆, γ 4 = 0). We expect that for small k z far away from the nexus there exists drumhead surface states bounded by the projected Dirac lines, and the number of these surface states for each k x and k z is determined by W (k x , k z ). This expectation is confirmed by the surface state spectrum shown in Fig. 10(a). On the other hand, we find that close to the nexus the surface states hybridize with bulk states similarly as in the case γ 3 > 0 [ Fig. 10(b),(c)]. On the other side of the nexus [Γ(k z ) < ∆/(2γ 1 )] one of the surface bands connects two electron-like bands to each other, and the other surface band connects two hole-like bands, so that the surface state spectrum is practically indistinguishable from the surface state spectrum in the case γ 3 > 0 [cf. Fig. 10(d) and Fig. 7