Appearance of flat surface bands in three-dimensional topological insulators in a ferromagnetic exchange field

We study the properties of the surface states in three-dimensional topological insulators in the presence of a ferromagnetic exchange field. We demonstrate that for layered materials like Bi$_2$Se$_3$ the surface states on the top surface behave qualitatively different than the surface states at the side surfaces. We show that the group velocity of the surface states can be tuned by the direction and strength of the exchange field. If the exchange field becomes larger than the bulk gap of the material, a phase transition into a topologically nontrivial semimetallic state occurs. In particular, the material becomes a Weyl semimetal, if the exchange field possesses a non-zero component perpendicular to the layers. Associated with the Weyl semimetallic state we show that Fermi arcs appear at the surface. Under certain circumstances either one-dimensional or even two-dimensional surface flat bands can appear. We show that the appearence of these flat bands is related to chiral symmetries of the system and can be understood in terms of topological winding numbers. In contrast to previous systems that have been suggested to possess surface flat bands, the present system has a much larger energy scale, allowing the observation of surface flat bands at room temperature. The flat bands are tunable in the sense that they can be turned on or off by rotation of the ferromagnetic exchange field. Our findings are supported by both numerical results on a finite system as well as approximate analytical results.


Introduction
A topological insulator is a material with an insulating energy gap in its bulk, but possesses conducting surface states due to significant spin-orbit coupling. The existence of these surface states is guaranteed by a topological invariant making them particularly robust against time-reversal invariant perturbations or disorder. Due to spin-orbit coupling the surface states are spin-momentum coupled allowing for interesting potential spintronics applications. The dispersion of the surface states forms a Dirac cone, i.e. the conduction electrons at the surface are effectively massless. This peculiar state of matter has first been suggested theoretically [1,2] and afterwards confirmed experimentally [3,4,5,6,7,8]. The number of materials identified as three dimensional topological insulators (3DTI) is steadily increasing [5,6,7,8,9,10,11,12].
In the present work we study 3DTI in the presence of a ferromagnetic exchange field. Experimentally, it has been demonstrated that such fields can be introduced into topological insulators either by doping with ferromagnetic dopants [15,16,17] or by proximity to ferromagnetic materials [18]. As the ferromagnetic exchange field breaks time-reversal symmetry, it allows for controlled modification or removal of the surface states and could lead to interesting effects or devices [13,14,15,19]. In previous work it has been pointed out that depending on the relative orientation of the exchange field with respect to the surface of the topological insulator, the surface states may open an energy gap or remain intact [20,21,22,23]. However, only small exchange splittings have been considered. In the present work we are considering exchange splittings up to the order of the gap of the topological insulator, which is up to 0.3 eV for present materials. Exchange splittings of the order of 1 eV can be reached with ferromagnetic materials.
In a recent work we have studied a two dimensional thin strip of a particle-hole symmetric topological insulator and found that under exchange fields of such strength an edge state flat band can appear [24]. Such flat bands are of particular interest, because the group velocity vanishes allowing highly localized wave packets. Flat bands have been found previously in other condensed matter systems like graphene, superfluid 3 He, or unconventional superconductors [25,26,27,28,29,30,31,32,33,34,35,36,37]. In particular, the appearence of flat bands in d-wave superconductors as surface Andreev bound states has been studied well in the past both theoretically and experimentally [33,34,35,38,39,40,41,42,43,44,45]. Such surface flat bands have been shown to lead to an enhanced barrier for vortex entry [46,47] or increased nonlinear electromagnetic reponse [48,49,50]. However, it remained an open question whether flat bands may also appear in three dimensional topological insulators under sufficiently strong ferromagnetic exchange fields.
In the present work, we will present a systematic investigation of the possible surface states of 3DTI and their behavior under a ferromagnetic exchange field. We will show under which circumstances surface flat bands appear. In particular, we identify a case in which a two dimensional flat band can be generated. We demonstrate that the appearence of our flat bands can be understood in terms of a classification recently proposed by Matsuura et al [51] using a topological invariant in the presence of a chiral symmetry. We will also show that the exchange field can produce highly anisotropic Dirac cones, i.e. that the group velocity is different in different directions. In this case the velocity can be tuned by rotation of the exchange field, i.e. rotation of the remanent magnetization of the ferromagnetism.
Recently the possibility of realizing a Weyl semimetallic state in pyrochlore iridates has raised a lot of interest [52]. This state is a generalization of the twodimensional Dirac electrons in graphene to a three-dimensional bulk system. In a Weyl semimetal conduction band and valence band touch each other only at a finite number of points. These so-called Weyl nodes are exceptionally stable for topological reasons. Of particular interest are the surface states of a Weyl semimetal which may form open Fermi "arcs" [52,53,54]. In the present work we show that a 3DTI with a sufficiently strong ferromagnetic exchange field becomes a Weyl semimetal in most cases. The surface flat bands are directly related to the appearence of surface Fermi arcs in this system. In contrast to previous proposals for surface flat bands in other systems like graphene, superfluid 3 He, or unconventional superconductors the present system has the advantage that the relevant energy scale is much larger (∼ 0.3 eV). This allows observation of the flat bands at room temperature, while for all other previous proposals cryogenic temperatures are necessary. In addition, the surface states can be tuned by rotation of the exchange field. For example, the flat band can be turned on and off by a rotation of the remanent magnetization of a ferromagnet. We demonstrate that such behavior can be achieved for realistic material parameters leading to new possible spintronic devices.

Models
As a starting point we consider the generic effective two-orbital Hamiltonian for a three dimensional topological insulator suggested already in several previous works [56,57,55,58]. To facilitate numerical calculations we choose the lattice regularized version that has been suggested by Li et al [55]: Here, ǫ 0 (k) = C + 2D 2 (1 − cos k x ) + 2D 2 (1 − cos k y ) + 2D 1 (1 − cos k z ), m 0 (k) = M − 2B 2 (1 − cos k x ) − 2B 2 (1 − cos k y ) − 2B 1 (1 − cos k z ), m 1 (k) = 2A 2 sin k x , m 2 (k) = 2A 2 sin k y , and m 3 (k) = 2A 1 sin k z . The Dirac Γ matrices are represented by Γ 0,1,2 = (I 2×2 ⊗ τ x , σ x ⊗ τ z , σ y ⊗ τ z ) in the spin-orbit basis. Here, the Pauli matrices in orbital space are denoted by τ i and the ones in spin space by σ i . As has been pointed out by Hao and Lee [59], there exist the following two different choices for Γ 3 : Γ 3 I = I 2×2 ⊗τ y and Γ 3 II = σ z ⊗τ z . These correspond to two different types of spin-orbit coupling in z-direction. We follow the convention of Hao and Lee and denote these two choices as model I and model II, respectively. Note that if k z = 0 there is no difference between the two models. Model I is isotropic within the xy-plane, but the coupling in z-direction is different. Thus one has qualitatively different behavior of surface states at a z-boundary than at an x-or y-boundary. For model II the spin-orbit coupling is isotropic and it is sufficient to consider surface states at one selected boundary, because the qualitative behavior is the same in all three spatial directions. It has been discussed in Ref. [57] that in the absence of an exchange field model I and model II are related by a unitary transformation and a 90 degree rotation within the xy-plane. However, this is not true anymore in the presence of an exchange field, because the spin operators are mapped to pseudospin operators under the unitary transformation as has been pointed out in Ref. [60]. Thus, for the purpose of the present work the two models become different in the presence of an exchange field. Model I is appropriate for Bi 2 Se 3 and its relatives.
The components of the ferromagnetic exchange field in x, y, and z-direction are denoted by V x,y,z , respectively, and are modeled by Zeeman terms in the Hamiltonian Eq. (1). The matrices for the exchange field components are given by Γ α = σ α ⊗ I 2×2 . The parameters A 1 , A 2 , B 1 , B 2 , C, D 1 , D 2 , and M have been derived from bandstructure calculations for the Bi 2 Se 3 family of materials in Refs. [56,57]. In our numerical calculations we will consider the case A i ≥ M > 0 and B i ≥ M as is relevant for these materials.

Symmetry considerations
Let us first discuss certain symmetries of the Hamiltonian (1) that will be of particular importance in this work. To begin with we focus on the particle-hole symmetric case C = D 1 = D 2 = 0, in which the effects become particularly clear and the topological invariant proposed by Matsuura et al [51] can be used. In Section 7 we will discuss the modifications that appear when particle-hole symmetry is slightly broken, as is the case in the Bi 2 Se 3 family of materials.
For C = D 1 = D 2 = 0 the four bulk bands for model I can be found by analytical diagonalization of the 4 × 4 matrices and are given by while for model II we have As is clear from these expressions, the bulk spectrum is fully symmetric around energy E = 0 for both models.
The Γ matrices introduced in the previous section respect the following commutation and anti-commutation relations where δ ij is Kronecker's delta symbol.
In the absence of an exchange field time reversal symmetry is respected. However, application of an exchange field in any direction breaks time reversal symmetry. According to the classification by Schnyder et al [61,62] a system with broken time reversal symmetry may still be topologically nontrivial, if a chiral symmetry is present. Let Θ be a chiral symmetry operator, which by definition anticommutes with the Hamiltonian, i.e.
Alternatively, we may also consider the symmetry operator Θ 3 = σ z ⊗ τ z , which happens to be identical to the operator Γ 3 II . Θ 3 anticommutes with Γ 0 , Γ 1 , Γ 2 , Γ 3 I , Γ x , and Γ y , but commutes with Γ 3 II and Γ z . Thus, for model I with exchange field within the x-y-plane Hamiltonian (1) with C = D 1 = D 2 = 0 possesses the chiral symmetry Θ 3 . For model II this symmetry is respected for k z = 0.
From these considerations we see that Hamiltonian (1) under certain circumstances falls into the AIII chiral symmetry class and we have identified three important symmetries.

Nonequivalent surface boundaries
Our aim is to calculate the energy dispersion of the surface states for both models, for all possible nonequivalent surfaces perpendicular to x-, y-, and z-directions, and for the corresponding directions of the exchange field. In this way we will determine all possible types of surface states that can appear in a 3DTI in a ferromagnetic exchange field.
For all cases we will present numerical calculations based on an exact diagonalization of Hamiltonian (1) on a finite size lattice of dimension 500 × 500 × 200. Periodical boundary conditions are employed parallel to the surface with 500 k-modes in both directions, while open boundary conditions are used perpendicular to the surface on 200 real space points. Our numerical results are compared with approximate analytical results for a continuous half space using a small k expansion of Hamiltonian (1) near the Γ point k = 0. For the appearance of the flat bands we will check our results using the topological winding number proposed by Matsuura et al [51].
In total we find that we need to consider seven nonequivalent cases: for model II the spin-orbit coupling in z-direction is of the same type as in x-and y-direction. For that reason it is sufficient to study a single boundary direction, which we choose to be a y-boundary, i.e. a boundary with y = const. As regards the direction of the exchange field we have to distinguish two nonequivalent cases here: parallel and perpendicular to the surface, i.e. V x = 0 and V y = 0. In contrast, for model I we have five nonequivalent cases. For model I the spin-orbit coupling in z-direction is of a different type than in x-and y-direction, but the in-plane coupling is still isotropic. Therefore we need to distinguish a z-boundary and a y-boundary. For the z-boundary there are again two nonequivalent directions for the exchange field: V x = 0 and V z = 0. For the y-boundary, however, all field directions are nonequivalent and we have three cases here.
We will see that among these seven cases there are three in which flat surface bands appear. One dimensional flat bands are found for model II with y-boundary and exchange field in x-direction as well as for model I with y-boundary and exchange field in z-direction. A two dimensional flat band is found for model I with y-boundary and exchange field in x-direction.

Model I
In this section we discuss the five nonequivalent cases for the particle-hole symmetric model I. We start with the more interesting case of a y-boundary.
5.1. Boundary perpendicular to the y-direction with finite V y In this case with V x = V z = 0 the bulk energy bands Eq. (2) simplify to the following expression: In the absence of an exchange field this band structure usually possesses a gap, because m 0 , m 1 , m 2 , and m 3 do not become zero simultaneously. Thus, the system is insulating. However, the gap closes when V y reaches a critical value V cr , which is derived in appendix I and is of the order of M . The Fermi surface at zero energy is then defined by the two equations m 1 = 0 and V 2 y = m 2 0 + m 2 2 + m 2 3 . These two equations define a line in three dimensional (k x , k y , k z ) space. Therefore the Fermi surface is one-dimensional and the system has entered a semimetallic state. If one looks at the Γ point k x = k y = k z = 0, where m 1 = m 2 = m 3 = 0 and m 0 = M it is clear that the semimetallic state is entered at V y = M or closely below. Thus, the parameter M > 0 sets the scale for the exchange field, at least for the range of the parameters A i and B i considered here. In appendix I we derive the ranges of the exchange field under which the system becomes semimetallic.
To find approximate analytical solutions for the surface states we expand Hamiltonian (1) up to second order in k y . If we assume a boundary in y-direction the momentum k y has to be replaced by the momentum operator −i∂ y . To find the surface states we search for nontrivial solutions of the Schrödinger equation that vanish both at y = 0 and for y → ∞. In this case the Hamiltonian can be written as where Here,m 0 (k) = M −2B 2 (1−cos k x )−2B 1 (1−cos k z ) and k = (k x , k z ). Γ 3 I anticommutes with H 0 . Γ y commutes with both H 0 and Γ 3 I , and Γ 1 anticommmutes with H 0 , Γ 3 I , and Γ y . In this case the eigenstates of H(k) are linear combinations of (up to four) eigenstates of H 0 (k) (see appendix II for a more detailed explanation). Surface states of H are superpositions of surface states of H 0 . Following the general procedure from appendix II we will first determine the surface states of H 0 and then deduce the ones of H from them.
For a system of finite width in y-direction, the energy of the surface states of H 0 behaves like e −L as a function of the system size L in the y-direction. Let us consider a half infinite system with a single boundary at y = 0. In this case the surface states of H 0 have zero energy. To find these zero energy eigenstates we can exploit that H 0 commutes with the operator Θ 4 = σ x ⊗ τ x . Then there exist simultaneous eigenstates of H 0 and Θ 4 . The eigenstates of Θ 4 are (1, 0, 0, 1) T and (0, 1, 1, 0) T with eigenvalue +1 and (1, 0, 0, −1) T and (0, 1, −1, 0) T with eigenvalue −1. We try the following two ansätze: where C is a normalization constant and f k (y) is solution of the equation The other two eigenstates of Θ 4 lead to exponentially increasing functions with y and thus cannot fulfil the boundary condition for y → ∞. Solving the differential equation (15) we find that f k (y) is given by This solution can only fulfil the boundary condition for y → ∞, ifm 0 (k) > 0. For those k values where this condition is not fulfilled anymore, a surface state does not exist.
Having determined the surface states of H 0 we can now infer the ones of H by noting that H ′ just couples the two solutions Eq. (13) and (14) as Here, u and v are the components of the solution spinor ξ k = (u k , v k ) T of the following eigenequation The full surface state solutions are then given by where with the corresponding eigenenergies These solutions show that we have surface states as long asm 0 (k) > 0. For small momenta the dispersion Eq. (20) shows that the presence of the exchange field in y-direction shifts the surface Dirac cone in k z -direction. The Dirac cone remains ungapped and its velocity is unchanged.
In Fig. 1 we show results obtained from numerical calculation of the eigenenergies of a finite slab on a finite size lattice for the present case showing nice agreement with the analytical result. Here and in the following we are showing results for the parameter choice A 1 = A 2 = B 1 = B 2 = M = 1 and C = D 1 = D 2 = 0. Results for parameters appropriate for Bi 2 Se 3 are discussed in section 7. Fig. 1(a) (for k z = 0) and (c) (for k x = 0) show the dispersion in the insulating state for a small exchange field of V y /M = 0.2. In Fig. 1 two dispersions are localized on the opposite surface, which is present in the numerical calculation. These states are related by parity to the ones found analytically above. Their dispersion is thus obtained from Eq. (20) by changing (k x , k z ) → (−k x , −k z ), i.e.
We note that the surface states Eq. (19) possess an interesting nontrivial spin texture. To see this we evaluate the expectation value of the spin components in the two orbitals. For orbital 1 the spin matrices can be written in the form where σ i for i ∈ {x, y, z} are the spin Pauli matrices. For orbital 2 we have analogouslŷ Using Eq. (19) we find the following expectation values in orbital 1: and in orbital 2: From these expressions we see that the spin rotates within the x-y plane. The spin direction of the two surface states is always opposite. The spin-x-component is opposite in the two orbitals, while the spin-y-component is the same in the two orbitals. Therefore, the total spin points in y-direction, perpendicular to the surface: Note, that while the total spin is perpendicular to the surface momentum, the partial spins in the two orbitals are not.

Boundary perpendicular to the y-direction with finite V x
Let us consider next the case that both V x and V y are nonzero, but still V z = 0. It is useful to go over to polar coordinates in this case and write V x = V 0 cos ϑ and V y = V 0 sin ϑ. The bulk energy bands Eq. (2) can then be brought into the following form: Compared with Eq. (9) this corresponds to a rotation of the exchange field within the x-y-plane. Again, the system is insulating in the absence of an exchange field and the gap closes, when V 0 reaches the critical value V cr . In the semimetallic state the Fermi surface is defined by the two equations m 1 sin ϑ − m 2 cos ϑ = 0 and V 2 0 = m 2 0 + m 2 3 + (m 1 cos ϑ + m 2 sin ϑ) 2 . Therefore, the Fermi surface is still onedimensional.
Determination of the surface states becomes more difficult now, because Γ x neither commutes nor anticommutes with H 0 in Eq. (11). As a result, it affects the spatial part of the surface states. We can, however, determine the surface states of the Hamiltonian H ′ 0 (k) still commutes with Θ 4 = σ x ⊗ τ x and we can thus look for zero energy states of H ′ 0 (k) using the same ansatz Eq. (13) and (14) as before. The surface state solutions of H ′ 0 (with a boundary at y = 0 and a half infinite system as before) are then found to be where C and C ′ are normalization constants. It is clear that state 1 exists only if m 0 (k) + V x > 0 and state 2 exists only ifm 0 (k) − V x > 0. The matrices Γ 1 , Γ 3 I , or Γ y neither commute nor anticommute with H ′ 0 . Thus they affect the spatial part of the surface states as well as the spin part. In this case we cannot separate the Hamiltonian into parts. However, if |k| and V y are small we can treat H ′ Eq. (12) as a perturbation using degenerate perturbation theory. We assume that the perturbation only couples the two surface states 1 and 2 to each other and neglect coupling to bulk states. This is a reasonable assumption, as the surface states are well localized and in most cases well separated in energy from the bulk states. In this case the overlap between the bulk states and the surface states is very small. Also, due to the symmetric energy spectrum around E = 0 for each bulk state with energy E there exists another one with energy −E. Thus, their contributions tend to cancel each other in perturbation theory.
If both states 1 and 2 exist, the surface states for the full Hamiltonian H are then approximately given by where The surface state eigenenergies in this case are found to be where is the spatial overlap of the two states Eq. (26) and (27). The dependence of β as a function of V x at k = 0 is shown in Fig. 2.
For small momenta the dispersion Eq. (29) shows that the presence of an exchange field within the xy-plane does not affect the presence of the surface Dirac cone. The Dirac cone is just shifted in k z -direction and remains ungapped. However, the velocity of the Dirac cone is isotropically suppressed by the factor β(V x , k). Thus, the presence of an exchange field component in x-direction allows tuning of the group velocity of the surface states, as shown in Fig. 2.
The spin texture of the surface states Eq. (28) turns out to be the same as in the previous section, except for the fact that all spin components are suppressed by the factor β(V x , k). Thus, the presence of the x-component of the exchange field V x leads to a suppression of the spin polarization of the surface states. Figure 3 shows dispersions with exchange field (V x , 0, 0). One can see from the figure that in this case k x -and k z -direction are equivalent. This is due to the fact that Γ x commutes with both Γ 1 and Γ 3 . In Fig. 1 (b) and (d) we see that a flat band appears, if V x > M . This flat band is apparently two dimensional, as it stays flat in both k x and k z direction. It still exists if we add a finite V y . The existence of this flat band goes beyond the perturbative treatment above. In the next subsection we discuss the existence of this flat band in terms of a topological invariant recently proposed by Matsuura et al. [51]

Existence of a two-dimensional flat band
Matsuura et al. [51] presented a general classification of the gapless topological phases like in semimetals or nodal superconductors. They showed that a generalized bulkboundary correspondence exists that relates the topological properties of the Fermi surface to the presence of protected flat bands at the surface of the system. In particular, it was found that the dimension of the surface flat band is always given by the dimension of the Fermi surface plus 1, if it exists (see Table V in Ref. [51]).
In the present case the system is topologically nontrivial and belongs to class AIII as discussed in section 3. The Fermi surface is one-dimensional and we may thus expect the appearance of a two-dimensional flat band.
The presence or absence of the flat band can be classified by a topological winding number. To construct this winding number, we use the chiral symmetry Θ 3 = σ z ⊗ τ z that was discussed in section 3 and is valid in the present case. Whenever a chiral symmetry is present, the Hamiltonian anticommutes with the symmetry operator. In this case the bulk Hamiltonian can be brought into off-diagonal block form by transforming to the eigenbasis of Θ 3 : where the block D(k) is found to be From the block D(k) we can define a winding number [64,51] Here, k ⊥ is the momentum component perpendicular to the surface. The winding number w is always an integer and depends on the momentum components parallel to the surface. It measures the phase change of the complex number det D(k), when k ⊥ runs through the Brillouin zone. If w(k || ) = 0, there exists no zero energy state for the momentum k || at the surface. If w(k || ) is nonzero, a zero energy surface state exists. In the present case, where we consider a boundary in y-direction we have Fig. 4 shows the phase of det D(k) for k x = 0 as a function of k y and k z in a color coded scale for the same set of parameters as in Fig. 3

(b) and (d).
Blue color corresponds to phase 0 and white to phase ±π. From the figure one notices that there exist a momentum space vortex and an anti-vortex at the positions (k y , k z ) = (0, ±1.318), around which the phase winds by 2π. These positions are points on the bulk Fermi surface. Note that det D(k) = 0 on the Fermi surface and the phase becomes singular there. From Fig. 4 it becomes clear that the winding number (33) becomes 1, when k z ∈ [−1.318, 1.318] and 0 outside. This is just the momentum range of the flat band seen in Fig. 3 (b). This example demonstrates that the winding number (33) correctly predicts the presence of the flat band.
Using Eq. (33) we can derive an analytical condition for the existence of the flat band, which is given in appendix III. In Figures 5 and 6 we show the regions in (k x , k z )-space in which the two-dimensional flat band on the y = 0 surface appears for different sets of parameters. These areas were calculated using the analytical condition from appendix III. Note, that the boundaries of these areas are just the projections of the one-dimensional Fermi surface onto the (k x , k z )-plane, as shown in appendix III. Figure 5 illustrates the effect of a rotation of the exchange field within the xy-plane on the flat band area. The flat band area is shown for exchange fields and V y = V 0 sin ϑ with V 0 = 2.8M and four angles of rotation ϑ. When the direction of the field is rotated from the x-direction into y-direction the size of the flat band area is reduced in x-direction but remains unchanged in z-direction. When the exchange field points in y-direction, the flat band finally disappears.
In appendix I we show that for M < 2A 2 1 /B 1 the minimal strength of the exchange field to create a semimetallic state and thus a two-dimensional surface flat band for the present case is given by V cr = M . However, if M > 2A 2 1 /B 1 (implicitly assuming B ≥ A and B ≥ M ), the minimal exchange field strength is given by the more complicated expression which is smaller than M . In this case we can have the situation that the flat band area is not simply connected anymore as shown in figures 6 (a) and (b). This case appears, when Eq. (112) from appendix I possesses two solutions instead of just one.

Boundary perpendicular to the y-direction with finite V z
Next, we consider the case that V z is nonzero and V x = 0. The component V y will be treated perturbatively like in section 5.2. The bulk energy bands Eq. (2) in this case simplify to Again, the system is insulating in the absence of an exchange field and the gap closes, when V z reaches the critical value V cr . In the semimetallic state the Fermi surface is defined now by three instead of two equations, m 1 = m 2 = 0 and V 2 z = m 2 0 + m 2 3 . For this reason, the Fermi surface becomes zero-dimensional, i.e. there are point nodes.
In this case the chiral symmetry Θ 3 does not hold anymore, because Θ 3 commutes with Γ z . However, as discussed in section 3 for the special case k x = 0 the system possesses the chiral symmetry Θ 1 . For that reason we may expect a one-dimensional flat band with k x = 0 in this case.
To determine the surface states in second order in k y one first notices that Γ z neither commutes nor anticommutes with H 0 . Therefore, it affects the spatial part of the surface states. However, similarly as in section 5.2, we can determine the zero energy surface states of the Hamiltonian This Hamiltonian commutes with the symmetry operator Θ 5 = σ z ⊗ τ x and it anticommutes with Θ 1 . As both Θ 1 and Θ 5 commute with each other, it is useful to look for surface state solutions among the common eigenstates of Θ 1 and Θ 5 .
These eigenstates are (−1, 1, 1, 1) T , (1, −1, 1, 1) T , (1, 1, −1, 1) T , and (1, 1, 1, −1) T . We thus try the following two ansätze: (the other two eigenstates lead to exponentially increasing functions again). We find that f k (y) is solution of the equations where the plus sign holds for ψ 1,k and the minus sign for ψ 2,k . Solving the differential equation (39) we find the solutions It is clear that state 1 exists only ifm 0 (k) + V z > 0 and state 2 exists only if Having determined the zero energy surface states of H ′′ 0 we can now try to obtain the ones of H = H ′′ 0 +m 1 Γ 1 +m 3 Γ 3 I +V y Γ y from them. First, one notices that Γ 1 = Θ 1 anticommutes with H ′′ 0 and the states (40) and (41) are already eigenstates of Γ 1 . Thus, these eigenstates are also eigenstates of H ′′ 0 + m 1 Γ 1 with energies ∓2A 2 sin k x . The operators Γ 3 I and Γ y neither commute nor anticommute with H ′′ 0 +m 1 Γ 1 . However, if k z and V y are small we can treat the terms m 3 Γ 3 I + V y Γ y perturbatively, again. For the same reasons as in section 5.2, we may assume that the perturbation only couples the two surface states 1 and 2 to each other. The surface states for the full Hamiltonian are then found to be of the form where The energies are given by Here, the spatial overlap β(V z , k) has the same functional form as in Eq. (30). For small momenta this dispersion shows that the y-component of the exchange field shifts the surface Dirac cone in k z -direction. The Dirac cone remains ungapped by the exchange field. The velocity of the Dirac cone is suppressed by the z-component of the exchange field only in k z -direction, but not in k x -direction. Thus, in the present case the exchange field can tune the group velocity of the surface electrons in an anisotropical way.
To determine the spin texture of the surface states we evaluate the spin expectation values in the two orbitals. For orbital 1 we find and in orbital 2: As in the previous cases the spin rotates within the x-y plane. The spin direction of the two surface states is always opposite. The spin-x-component is opposite in the two orbitals, while the spin-y-component is identical. With increasing V z the spin-ycomponent is suppressed by the β(V z , k) factor. The total spin points in y-direction again, perpendicular to the surface: Figure 7 shows results for the energy dispersions with finite V z obtained from numerical calculations on a finite size system. One can see from the figure that the directions k x and k z differ. This is again due to the fact that Γ z commutes with Γ 3 I but anticommutes with Γ 1 . We see from the figure that a one dimensional flat band

Existence of a one-dimensional flat band
Similarly as in the case with nonzero V x the appearance of this one-dimensional flat band can be understood using the classification of Matsuura et al. [51] For that purpose we consider the Hamiltonian without the Γ 1 term: If one considers this Hamiltonian for k x = 0 as a function of the two coordinates k y and k z , its Fermi surface will be a point node in the two-dimensional Brillouin zone. However, one can also consider this Hamiltonian as a function of the three coordinates k x , k y , and k z , keeping the k x dependence in m 0 (k). Then, its Fermi surface will be a line node in the three-dimensional Brillouin zone. In any case, H 1 belongs to class AIII due to the chiral symmetry Θ 1 . Using the chiral symmetry Θ 1 we can bring H 1 into off-diagonal block form, similarly as in section 5.3: where the block D 1 (k) is found to be Using this block, we can again define the winding number Following the method outlined in appendix III we can derive an analytical condition for the existence of a flat band of H 1 which reads For k x = 0 this yields a range of k z values for which a one-dimensional flat band exists, consistent with the numerical result in Fig. 7. If we consider H 1 as a function of three coordinates k x , k y , and k z , this condition tells us that H 1 actually possesses a two-dimensional surface flat band within a certain area in (k x , k z )-space. Now, H = H 1 + m 1 Γ 1 and Γ 1 anticommutes with H 1 . This means that the zero energy states of H 1 are eigenstates of Γ 1 = Θ 1 , too. As a result, the dispersion of the surface flat band of the full Hamiltonian H is given by Thus, the surface state dispersions become highly anisotropic, being flat in k zdirection, but dispersive in k x -direction. In Figure 8 we illustrate an example of the general case, where all three components of the exchange field are nonzero. One sees from the figure that at small fields the surface states are splitted, similarly as in Figure 1 (a) and (c). However, when the magnitude of the exchange field exceeds M , a one-dimensional flat band appears in k z direction, similarly as in Figure 7 (b) and (d).
In order to understand the appearance of a flat band in this general case let us consider the following partial bulk Hamiltonian: We first construct a symmetry operator Θ 13 that anticommutes with this H 0 . We first note that the operator Θ 1 anticommutes with Γ 0 , Γ 2 , Γ y , and Γ z , but commutes with Γ x , while the operator Θ 3 anticommutes with Γ 0 , Γ 2 , Γ x , and Γ y , but commutes with Γ z . If we choose, however, the following linear superposition of Θ 1 and Θ 3 it is easy to see that Θ 13 anticommutes with the four operators Γ 0 , Γ 2 , Γ y , and Thus, the partial Hamiltonian H 0 possesses a chiral symmetry Θ 13 , which depends on the direction of the exchange field. For V z = 0 the symmetry operator reduces to Θ 3 , corresponding to the case in section 5.3 and for V x = 0 the symmetry reduces to the case discussed in section 5.5.
Using the chiral symmetry Θ 13 we can bring H 0 into off-diagonal block form by transforming to the eigenbasis of Θ 13 . There are two eigenvectors of Θ 13 with eigenvalue -1, which are η 1 = (0, 0, sin χ 2 , cos χ 2 ) T and η 2 = (cos χ 2 , − sin χ 2 , 0, 0) T and two eigenvectors with eigenvalue +1, which are η 3 = (sin χ 2 , cos χ 2 , 0, 0) T and η 4 = (0, 0, cos χ 2 , − sin χ 2 ) T . In this basis H 0 becomes block off-diagonal where the block D 2 (k) is and its determinant The Hamiltonian H 0 possesses a two-dimensional zero energy surface flat band, if the corresponding winding number becomes nonzero. Following appendix III we find the criterion that In the vicinity of k x = k z = 0 this criterion can usually be fulfilled for V > M , if V does not become too large.
Having seen that the partial Hamiltonian H 0 possesses a two-dimensional zero energy surface flat band we can apply as a perturbation to find the approximate surface state dispersion for the full Hamiltonian. For that purpose we need the surface state wave function ψ s . As the zero energy state is a common eigenstate of H 0 and Θ 13 , ψ s will be a superposition of two eigenvectors of Θ 13 with the same eigenvalue, like The two functions f 1 (y) and f 2 (y) have to be determined from the differential equation H 0 (k y → −i∂ y )ψ s = 0 and are found to be superpositions of three exponentially decaying functions. An analogous superposition of the +1 eigenstates of Θ 13 leads to exponentially increasing functions, which cannot fulfil the boundary conditions. They correspond to solutions localized at the opposite boundary. Both functions f 1 (y) and f 2 (y) fulfil f 1 (y = 0) = f 2 (y = 0) = 0 and are normalized such that The energy of the perturbed system is given by By direct calculation we find that and As a result we find for the energy of the surface states of the full Hamiltonian: From this expression we see that for V z = 0 we always have a one-dimensional flat band in k z -direction. The group velocity in x-direction can be tuned by rotating the exchange field within the xz-plane and vanishes when V z becomes zero. This expression shows how the one-dimensional flat band develops into the two-dimensional flat band for exchange field within the xy-plane.

Weyl semimetal
In this section we show that the semimetallic state of model I for an exchange field V > V cr and V z = 0 is actually a realization of a Weyl semimetal. The Weyl semimetallic phase can be viewed as a three-dimensional generalization of the twodimensional Dirac electrons in graphene, as has been pointed out recently [52,53,54]. In contrast to graphene, just two linearly dispersing adjacent bands touch at a finite number of points in the three-dimensional Brillouin zone. In the vicinity of these Weyl nodes, which have also been termed "diabolic" points [66], the effective twoband Hamiltonian can be written in the form where the Pauli matrices σ i do not need to refer to the spin degree of freedom. Such a diabolic point is exceptionally stable due to topology, as arbitrary perturbations cannot remove it unless an annihilation with another diabolic point occurs. A recent theoretical work proposed the appearence of a Weyl semimetallic phase in pyrochlore iridates [52]. It was shown that the surface states of the Weyl semimetal may form open Fermi "arcs", i.e. Fermi lines which terminate at the projection of the diabolic points onto the surface Brillouin zone.
In appendix I we showed that for V > V cr model I enters a semimetallic phase. If V z = 0 the bulk spectrum indeed possesses either two or four Fermi points on the k z -axis, where two bands touch each other. We still need to show that around these points the Hamiltonian can be written in a form like Eq. (62). In order to do so we use a similar technique as we have used for determination of the dispersion of the surface states. Let k 0 = (0, 0, k z,0 ) be the position of a Fermi point. We first consider the partial bulk Hamiltonian at k 0 and construct a symmetry operator Θ 12 that anticommutes with it. The zero energy eigenstates of H 0 are then simultaneous eigenstates of Θ 12 . From the eigenstates of Θ 12 we determine the two zero energy eigenstates of H 0 at the Fermi point. We then expand the full Hamiltonian to lowest order in k x , k y , and k z − k z,0 around the Fermi points. This leads to a perturbation of the form The effective low energy 2 × 2 Hamiltonian near the Fermi points is obtained by degenerate perturbation theory within the subspace of the two said eigenstates.
To construct the symmetry operator Θ 12 we first note that the operator Θ 1 anticommutes with Γ 0 , Γ 3 I , Γ y , and Γ z , but commutes with Γ x , while the operator Θ 2 anticommutes with Γ 0 , Γ 3 I , Γ x , and Γ z , but commutes with Γ y . If we choose the following linear superposition of Θ 1 and Θ 2 it is easy to see that Θ 12 anticommutes with the four operators Γ 0 , Γ 3 I , Γ z , and V x Γ x + V y Γ y . There are two eigenvectors of Θ 12 with eigenvalue -1, which are η 1 = (0, 0, ie −iϕ , 1) T and η 2 = (−ie −iϕ , 1, 0, 0) T and two eigenvectors with eigenvalue +1, which are η 3 = (0, 0, −ie −iϕ , 1) T and η 4 = (ie −iϕ , 1, 0, 0) T . One of the two zero energy eigenstates of H 0 is a linear combination of η 1 and η 2 , while the other one is a linear combination of η 3 and η 4 . After a staightforward calculation which exploits the fact that at the Fermi points we have m 2 0 + m 2 3 = V 2 x + V 2 y + V 2 z (see appendix I), we find the following eigenstates of H 0 : where e iχ = m 0 + im 3 To find the effective low energy 2×2 Hamiltonian near the Fermi points for convenience we use the two basis states In this basis the Hamiltonian H ′ becomes Apparently, this is an anisotropic Weyl-type Hamiltonian. For example for V x = 0 this expression simplifies to Having seen that model I for V > V cr and V z = 0 is a Weyl semimetal we can see now that the one-dimensional surface flat band from the previous section Eq. (49) is just a Fermi "arc" in the sense of Ref. [52]. It exists only in a limited range of k z values given by Eq. (48). The end points of the Fermi arc are just the projections of the Weyl nodes onto the surface Brillouin zone (k x , k z ), as one sees by setting k x = 0 in Eq. (48) and comparison with Eqs. (108) and (111) in appendix I.

Boundary perpendicular to the z-direction
The case with a z-boundary is much easier to treat than the case with a y-boundary. In this case we can decompose the Hamiltonian in the following way: Here,m 0 (k) = M − 2B 2 (1 − cos k x ) − 2B 2 (1 − cos k y ) and k = (k x , k y ) now.
We are now in the lucky situation that Γ x , Γ y and Γ z all commute with H 0 , and both Γ 1 and Γ 2 anticommute with it. Thus, the surface states for an exchange field in arbitrary direction can be found from the zero energy surface states of H 0 using the method from appendix II.
We first determine the zero energy surface states of H 0 by noting that H 0 commutes with Γ z and anticommutes with the operator I 2×2 ⊗ τ z . The common eigenstates of these two operators are just (1, 0, 0, 0) T , (0, 1, 0, 0) T , (0, 0, 1, 0) T , and (0, 0, 0, 1) T . We try the following two ansätze (the other two eigenstates leading to exponentially increasing functions again): where f k (z) is solution of the equation Solving the differential equation (76) we find that f k (z) is given by This solution can only fulfil the boundary condition for z → ∞, ifm 0 (k) > 0. For those k values where this condition is not fulfilled anymore, a surface state does not exist. The form of the solutions Eq. (74) and (75) means that the surface state only occupies orbital 1, leaving orbital 2 empty. Correspondingly, the surface state at the opposite surface of the system is found to occupy orbital 2 only. To obtain the surface states of H from the ones of H 0 we have to determine those linear combinations of ψ 1,k and ψ 2,k that diagonalize H ′ , i.e.
If we write ξ i (k) = (a i (k), b i (k)) T the coefficients are determined from the equation The eigenenergies are given by This expression tells us that the components of the exchange field parallel to the surface (V x and V y ) for small momenta just shift and split the surface Dirac cone without opening a gap and without changing the group velocity. The component V z perpendicular to the surface opens a gap, however. This behavior is in agreement with previous work [20]. A flat band does not appear in this geometry. Even though the system is still a Weyl semimetal in the bulk, a surface Fermi arc does not appear, because the Weyl nodes all sit on the k z -axis. Thus, their projection onto the surface Brillouin zone is the single point k x = k y = 0 and the Fermi arc is not present. The full surface state wave functions can be written in the form Here, we have introduced two spherical angles φ and θ that define the direction of the vector (m 1 + V x , m 2 + V y , V z ): For the spin texture of these surface states we find Ψ ±,k |ŝ 1,x | Ψ ±,k = ± cos φ k sin θ k Ψ ±,k |ŝ 1,y | Ψ ±,k = ± sin φ k sin θ k Ψ ±,k |ŝ 1,z | Ψ ±,k = ± cos θ k The spin in orbital 2 vanishes. This means that the spin is directed along the vector (m 1 + V x , m 2 + V y , V z ). For V z = 0 the spin-z component vanishes and the spin is oriented within the plane of the surface, in contrast to the cases with the y-boundary. Figure 9 shows numerical dispersions from a finite size system with finite V x in agreement with Eq. (81). For larger values of the exchange field V x > M the bulk gap has closed. The surface states then exist within the bulk bands. Figure 10 shows numerical dispersions with finite V z confirming the opening of a gap for small values of V z . Again, for larger values of the exchange field V z > M the surface states only exist within the bulk bands and no flat band appears.

Model II
In this section we discuss the two nonequivalent cases for the particle-hole symmetric model II. Model II is distinguished from model I only in the coupling in k z -direction by the matrix Γ 3 II = σ z ⊗ τ z instead of Γ 3 I . The other Γ-matrices are the same as in model I. The matrix Γ 3 II has different commutation and anti-commutation relations than Γ 3 I . Model II is more symmetric in the sense, that the three spatial directions possess equivalent couplings. Therefore, it does not matter which direction of the boundary we consider. For convenience we choose a boundary in y-direction, because this allows us to build on the results from model I found in the previous section.

Finite V x and V y
The case with both V x and V y nonzero, but V z = 0 can be treated along the same lines as has been discussed for model I in section 5.2. We can again go over to polar coordinates in this case and write V x = V 0 cos ϑ and V y = V 0 sin ϑ. The bulk energy bands Eq. (3) can then be brought into the form: Again, the system is insulating in the absence of an exchange field and the gap closes, when V 0 reaches a critical value V cr ∼ M . In contrast to model I, in model II V cr depends on the direction of the exchange field as detailed in appendix I. In the semimetallic state the Fermi surface is defined by three equations m 1 sin ϑ−m 2 cos ϑ = 0, m 3 = 0, and V 2 0 = m 2 0 + (m 1 cos ϑ + m 2 sin ϑ) 2 now. Therefore, the Fermi surface is pointlike, i.e. zero-dimensional. Similarly as in section 5.6 it can be shown by linearization around the point nodes that the system is a Weyl semimetal in this case, too.
To determine the surface states, we can start from the same Hamiltonian H ′ 0 in Eq. (25) and its zero energy surface states in Eq. (26) and (27). We can then treat the terms using degenerate perturbation theory again. In contrast to the matrix Γ 3 I the matrix Γ 3 II leads to a diagonal coupling of the two surface states instead of an off-diagonal coupling. As a result the dispersion of the surface states in k z -direction remains unaffected by the spatial overlap β(V x , k) Eq. (30). Consequently, we find the following surface state dispersions for the full Hamiltonian within perturbation theory: This dispersion shows that the V y component of the exchange field, which is perpendicular to the surface, opens a gap in the surface state spectrum for small values of V y . The component V x parallel to the surface leads to an anisotropic reduction of the dispersion in k x -direction, but not in k z -direction. In Figures 11 and 12 we show the corresponding numerical results on a finite lattice for V x nonzero and V y nonzero, respectively, which agree with this behavior. The surface state wave functions for this case can be written in the form Here, we have introduced two spherical angles φ k and θ k that define the direction of the vector (βm 1 , βV y , m 3 ): For the spin texture of these surface states we find in orbital 1: In orbital 2 the x-and z-components of the spin turn out to be inverse: Here, we see that in contrast to the corresponding case for model I the spin possesses components in all three spatial directions. In the limit V y → 0, the angle φ k goes to zero and the spin-y component vanishes. The spatial overlap factor β is seen to suppress only the x-and y-components of the spin. For the total spin we find: The total spin is directed perpendicular to the surface. It vanishes in the limit V y → 0. From the numerical results in Fig. 11 (b) and (d) we see that a one-dimensional flat band appears for nonzero V x > M . Again, the appearence of this flat band can be understood by a topological winding number. In the present case the system obeys the chiral symmetry Θ 3 for k z = 0, as was discussed in section 3. The Fermi surface is zero dimensional, so we may expect a one-dimensional flat band. As the Hamiltonians for model I and model II are identical for k z = 0, the off-diagonal block form of the Hamiltonian is the same as in Eq. (31), i.e.
The corresponding winding number is then given by The one-dimensional flat band in the present case thus corresponds to a k z = 0 cut of the flat band area shown in Fig. 5 (a). We can find the full dispersion of this flat band also for finite k z by noting that the matrix Γ 3 II anticommutes with the Hamiltonian for k z = 0. As a result the zero energy states for k z = 0 are eigenstates of Γ 3 II , too. Therefore, the dispersion of the surface flat band for finite k z is given by Like for model I the zero energy surface of this one-dimensional flat band is a Fermi arc, whose end points are the projections of the Weyl nodes onto the surface Brillouin zone.

Finite V z and V y
The case with both V z and V y nonzero, but V x = 0 can be treated like the corresponding case for model I in section 5.4. The Fermi surface turns out to be zero-dimensional and the system possesses the chiral symmetry Θ 1 for k x = 0. For determination of the surface states we can start from the zero energy surface states Eqs. (40) and (41) of the Hamiltonian Eq. (36) and treat H ′ = m 1 Γ 1 + m 3 Γ 3 II + V y Γ y as a perturbation. The energies of the surface states are then found to be This is the same kind of dispersion as in Eq. (86) with the roles of the x-and zcoordinates interchanged. The spin texture of the surface states in orbital 1 is found to be: and in orbital 2: For the total spin we find: The total spin is again directed perpendicular to the surface and vanishes in the limit V y → 0. We find a one-dimensional flat band for finite V z > M , that can be understood by a topological winding number using the chiral symmetry Θ 1 for k x = 0, analogously to the previous case. As the matrix Γ 1 anticommutes with the Hamiltonian for k x = 0 we eventually find the dispersion of the surface flat band in this case as Again, the zero energy surface of this one-dimensional flat band is a Fermi arc, whose end points are the projections of the Weyl nodes onto the surface Brillouin zone.

Effect of broken particle-hole symmetry
So far we have studied the case C = D 1 = D 2 = 0. When these parameters become nonzero, the particle-hole symmetry of the system is broken. Also, the chiral symmetries discussed in section 3 are not obeyed anymore. For this reason the topological winding numbers that we used in the previous sections to determine the presence or absence of a surface flat band cannot be used anymore. However, this does not mean that the surface bands completely disappear. In the following we will demonstrate by both numerical and analytical calculation that surface bands, which are energetically well separated from the bulk bands still exist in the broken particle-hole case. We will see that the surface bands become dispersive now with the dispersion increasing proportional to D 1 and D 2 . Similar behavior has been noted in other systems before as well [51,30,32]. For the numerical calculations we use parameters that are realistic for Bi 2 Se 3 and have been given in Ref. [56]. The Hamiltonian is now Here, we have used the lattice constants a = 4.14Å and 15c = 28.64Å for conversion into our lattice model. In the following, we use the same parameters for both model I and II. It is clear that the bulk energy bands are just shifted by ǫ 0 (k) with respect to the particle-hole symmetric cases studied in the previous sections. The former zero energy points, at which two bulk bands touch each other, will then have energy ǫ 0 (k). For model I with V z = 0 and V > V cr we had a one-dimensional Fermi surface, whose degeneracy will now be lifted. Thus, these nodes will generally overlap with the two touching bulk bands. However, for the particle-hole symmetric model I with V z = 0 and V > M we have a Weyl semimetal with just two Fermi points. As these two Weyl nodes sit at symmetry related momentum points, they are shifted by the same amount ǫ 0 (k). Therefore, the particle-hole broken system remains a Weyl semimetal unless the dispersion ǫ 0 (k) becomes so large that the Weyl nodes start to overlap with the bulk bands. This same argument also holds for model II with V > M . Thus, the Weyl semimetallic phases in model I and II are preserved under not too large particlehole symmetry breaking and we may still expect the presence of surface Fermi arcs, as we will show explicitly from our numerical calculations below. In the following we restrict the discussion to the geometries in which we found surface flat bands for the particle-hole symmetric cases.

Model I with boundary perpendicular to the y-direction
In the absence of an exchange field the surface state dispersions can be found by replacing the momentum k y by the momentum operator −i∂ y and solving the 4 × 4 matrix Schrödinger equation directly with an exponential ansatz following Ref. [58]. This way one finds where When an exchange field is turned on, the Schrödinger equation leads to 8th order polynomials, whose zeroes cannot be given in closed form. However, analytical results can be obtained for fields in high symmetry directions and small values of momentum k x and k z by expansion.
When the exchange field points into y-direction we find the following dispersions Here, β 3 is a spatial overlap factor of the form Eq. (96) tells us that the Dirac cone remains ungapped for an exchange field in ydirection and is shifted in k z direction, like in the particle-hole symmetric case in section 5.1. With exchange field in x-direction the following low field dispersions are found where the spatial overlap In contrast to Eq. (29) we now find the opening of a gap in the dispersion, which is of the order of t 2 V x . Corresponding numerical results on a finite lattice for V x = 0.2M are shown in Fig. 13 (a) and (c). If we increase V x beyond M , we enter a state that corresponds to the two-dimensional flat band state discussed in section 5.3. As one can see in Fig. 13 (b) and (d) the surface band becomes dispersive now. However, it still exists and remains well separated from the two bulk bands. The dispersion of the surface band over the surface Brillouin zone is shown in color coded scale in Fig. 14. The dispersion is much smaller along k x -direction than in k z -direction. An approximate analytical expression valid for small k x and k z can be obtained from Eq. (98), if one sets β 4 = 0. In the total density of states of this system the surface band appears as a peak, similar to what was found for the edge states in a twodimensional system in our previous work (see Fig. 5 in Ref. [24]).
With exchange field in z-direction we find the dispersions where β 5 is given by . From this dispersion we see that the Dirac cone remains ungapped for an exchange field in z-direction and is shifted in k x direction. The numerical results shown in Fig. 15 (a) and (c) confirm this behavior. In Fig. 15 (b) and (d) results are shown for V z > M . From these figures we see that also the one-dimensional flat band discussed in section 5.5 becomes dispersive now due to the broken particle-hole symmetry. However, the system remains a Weyl semimetal, as was pointed out above. The projections of the two Weyl nodes onto the surface Brillouin zone for V z = 4M are found at (k x , k z ) = (0, ±0.66) in Fig. 15 (d).
These are the points, where the surface band ends. In Fig. 16 we show the lines of constant energy E = 0.144 eV for the states on the y = 0 surface (solid line) and on the y = L surface (dashed line). Here, one can see that the line of constant energy at the energy of the Weyl nodes becomes a curved and open Fermi arc, as expected in a Weyl semimetal [52,53]. In contrast, the Fermi arc was a straight line in the particle-hole symmetric case discussed in section 5.6.

Model I with boundary perpendicular to the z-direction
In this case we find approximate surface state dispersions for small momenta k x and k y of the form where t 1 = D 1 /B 1 . Here, the z-component of the exchange field opens a gap in the Dirac dispersion, while both x-and y-components lead to a shift of the Dirac cone, leaving it intact. In this geometry there appear no Fermi arcs at the surface, as in the corresponding particle-hole symmetric case. This can be understood from the fact that both Weyl nodes sit on the k z -axis. For that reason their projections onto a k z = const. plane fall on top of each other and the Fermi arc shrinks to zero. It is interesting to see that the physics of the Weyl semimetal state cannot be observed on a z-surface due to the structure of the Weyl state here. Experimentally one needs to look at the side surfaces or at a surface having a finite angle with the z-surface to observe the Fermi arc.

Model II with boundary perpendicular to the z-direction
In the absence of an exchange field the surface state dispersions are given by where t 1 = D 1 /B 1 ≈ 0.13. If t 1 > 1 surface states do not exist. When the exchange field points into x-direction we find the following surface state dispersions where β 6 is given by . From this dispersion we see that the Dirac cone remains ungapped for an exchange field in x-direction and is shifted in k y direction. Corresponding numerical results are shown in Fig. 17 (a) and (c) confirming this behavior. In Fig. 17 (b) and (d) we show results for V x > M . Here, we see that again the corresponding one-dimensional flat band discussed in section 6.1 becomes dispersive now due to the broken particle-hole symmetry. Again, this system remains a Weyl semimetal, however. As pointed out in appendix I, in model II the Weyl nodes follow the direction of the exchange field and sit on the k x -axis in the present case. The projections of the two Weyl nodes onto the surface Brillouin zone for V z = 4M are found at (k x , k y ) = (±0.623, 0) in Fig. 17 (d). These are the points, where the surface band ends. Similarly as in Fig. 16 Thus, the application of an exchange field perpendicular to the surface leads to a gap in the Dirac dispersion. When V z is increased beyond V cr no flat band appears. This can be understood from the fact that in model II the Weyl nodes follow the field direction. Thus, for fields perpendicular to the surface their projections onto the surface Brillouin zone fall on top of each other and the Fermi arc disappears.

Summary and Conclusions
We have studied the modification of the surface states of a three dimensional topological insulator by a ferromagnetic exchange field using the two models that are presently discussed for topological insulators. For model I, which is appropriate for the Bi 2 Se 3 class of layered materials the surface states on a side surface behave qualitatively different than the ones on the top or bottom surface. For exchange fields smaller than the bulk gap the velocity of the Dirac cone can be tuned down to smaller values in an anisotropic way depending on the direction of the exchange field. For exchange fields larger than a critical value of the order of the bulk gap the system becomes a topologically nontrivial semimetal.
We have shown that in a particle-hole symmetric system the Fermi surface of this semimetal is a line in momentum space, if the exchange field is directed within the xy-plane. In this case a two-dimensional flat band appears at a side surface. If the exchange field possesses a finite component in z-direction, there exist only single points in the Brillouin zone, where the bulk gap vanishes. We have shown that in this general case the system becomes a Weyl semimetal. Associated with this peculiar state of matter we find Fermi arcs at the side surfaces and one-dimensional flat bands.
If particle-hole symmetry is not obeyed, the flat bands become dispersive, but remain well separated from the bulk bands. The Weyl semimetallic phase and the existence of a surface Fermi arc is preserved also under broken particle-hole symmetry.
In model II, which is more isotropic than model I, the behavior on the top and side surfaces is qualitatively the same. For small exchange fields, it is again possible to tune the velocity of the Dirac cone by the exchange field anisotropically. For exchange fields larger than a critical value of the order of the bulk gap model II always enters a Weyl semimetal state. Associated with this we also find Fermi arcs and one-dimensional surface flat bands in this case. However, in model II there is no case in which a two-dimensional surface flat band appears.
If particle-hole symmetry is not obeyed, the flat bands become again dispersive. The Weyl semimetallic phase and the existence of a surface Fermi arc again is preserved under broken particle-hole symmetry.
The Weyl nodes in model I sit on the k z axis, while they move with the direction of the exchange field in model II. As a consequence, in model I the Fermi arcs and flat bands do not appear on the top and bottom surfaces, while in model II this is possible.
We have shown that the appearence of flat bands in our particle-hole symmetric cases could always be classified by topological invariants that have recently been given by Matsuura et al [51]. The invariants are related to different chiral symmetries in the different cases.
Surface flat bands have been proposed in systems like graphene, superfluid 3 He, or unconventional superconductors before. However, in these systems cryogenic temperatures are required to observe the flat bands. In the materials discussed here, the energy scale is set by the bulk gap, which is of the order of 0.3 eV in Bi 2 Se 3 . Thus, flat bands could be observed already at room temperature. The fact that one can turn on or off a surface flat band by rotation of the magnetization of a ferromagnet makes the present system particularly interesting for applications in spintronics.
For model I we start from the bulk bandstructure of the four bands given in Eq. (2). We omit the diagonal term ǫ 0 (k)I 4×4 in the Hamiltonian as it just shifts the four bands by ǫ 0 (k) and thus does not affect the direct gap of the system. For large values of D 1 and D 2 this term can lead to an indirect gap, however.
As is clear from the symmetric spectrum Eq. (2) the gap closes when there exist momenta k for which E I i (k) becomes zero. To determine those momenta it is beneficial to write the exchange field in spherical coordinates, i.e. V x = V cos ϕ cos ϑ, V y = V sin ϕ cos ϑ, and V z = V sin ϑ. If one exploits the fact that can then only be fulfilled for selected values of k z . Thus the system will have single Fermi points, when a solution exists. For the special case ϑ = 0, i.e. when the exchange field lies within the x-y-plane, there are just two equations to be fulfilled, i.e. m 1 sin ϕ = m 2 cos ϕ and V 2 = m 2 0 + m 2 3 + (m 1 cos ϕ + m 2 sin ϕ) 2 . In this case the Fermi surface will be a line in momentum space.
To determine the ranges of the exchange field for which these solutions exist, we determine the minimum and maximum value of m 2 0 + m 2 3 . Let us set c = cos k z . Then we have to minimize the function in the interval c ∈ [−1, 1]. Here, we have set for k x = k y = 0 M − 4B 2 for (k x = 0 and k y = π) or (k x = π and k y = 0) M − 8B 2 for k x = k y = π The function f (c) is quadratic in c, so there can only be a single extremum within c ∈ [−1, 1] or the extremum will be on the boundaries c = ±1. On the boundaries we have As we assume B 1 > M > 0 and B 2 > M > 0, we have f (c = −1) > f (c = 1). To determine a possible extremum inside the interval c ∈ [−1, 1] we take the derivative of f (c) yielding The second derivative of f (c) is As we assume B 1 > A 1 , this is positive. We thus find that there exists a minimum inside the interval c ∈ [−1, 1], if which is equivalent to the condition This condition can only be fulfilled for positive M ′ , i.e. only for k x = k y = 0. Thus, in the case M B 1 ≥ 2A 2 1 , the minimum for the exchange field can be found by introducing Eq. (109) into Eq. (108). After some algebra one finds Thus we find that the minimum critical value V cr for the magnitude of the exchange field, which is necessary to bring the system into the semimetallic phase, is given by In total we find the following three ranges for the magnitude of the exchange field, in which the system becomes semimetallic: . For B 2 ≥ B 1 these ranges overlap, while for B 2 < B 1 they are separate. As B 1 and B 2 are usually of the order of one to several eV, we do not expect that large enough exchange fields can be applied in practice to actually observe these different ranges. However, the minimum critical field V cr , which is of the order of M or less, is within experimental reach. The position of the Fermi points is found from the quadratic equation For V ∈ [V cr , 4B 1 − M ] the Fermi points sit on the k z -axis (k x = k y = 0). Their positions are given by For the special case V z = 0 the minimum value of the exchange field to bring the system into the semimetallic state remains the same. In this case the function V 2 = m 2 0 + m 2 3 + (m 1 cos ϕ + m 2 sin ϕ) 2 has to be minimized. However, as the minimum of m 2 0 + m 2 3 occurs at k x = k y = 0 and m 1 = m 2 = 0 there, the minimum of V remains the same value V cr Eq. (110). Also, for V > V cr the two (or four) points determined above still lie on the Fermi surface. However, as pointed out above, the Fermi surface becomes a line (or two lines) now, which approximately run within the plane perpendicular to the exchange field (The precise condition is that the vector (m 1 , m 2 ) T = 2A 2 (sin k x , sin k y ) T should be perpendicular to the exchange field (V x , V y ).) Let us next look at model II. We again write the exchange field in spherical coordinates, i.e. V x = V cos ϕ cos ϑ, V y = V sin ϕ cos ϑ, and V z = V sin ϑ. This time Eq. (3) can be written in the form E II i (k) = ± (m 1 sin ϕ − m 2 cos ϕ) 2 + (m 1 cos ϕ sin ϑ + m 2 sin ϕ sin ϑ − m 3 cos ϑ) 2 + + V ± m 2 0 + (m 1 cos ϕ cos ϑ + m 2 sin ϕ cos ϑ + m 3 sin ϑ) This expression can become zero only, if all three squared expressions below the squareroot vanish simultaneously. This means that the system will have Fermi points, when a zero energy solution exists. The first two expression become zero, if the vector (m 1 , m 2 , m 3 ) is parallel to (cos ϕ cos ϑ, sin ϕ cos ϑ, sin ϑ), i.e. parallel to the exchange field. The third equation can then be written as It is clear that the minimum critical value V cr for the magnitude of the exchange field, which is necessary to bring the system into the semimetallic phase fulfils V cr ≤ M , because a zero energy state can always be found for k x = k y = k z = 0 and V = M . A general expression for V cr can in principle be obtained analytically for exchange field in any direction. However, the expressions become quite complicated. Therefore, here we focus on the two cases that the exchange field points either in x-direction or in z-direction.
For exchange field in x-direction we have m 2 = m 3 = 0. The minimum can thus be found on the k x -axis by minimizing the function f (c) Eq. (108) with c = cos k x and B 1 and A 1 being replaced by B 2 and A 2 . Following the same calculation as for model I we then find for the critical value V cr,x in this case: Analogously to model I, for V > V cr,x we have either two or four Fermi points lying on the k x -axis. Their positions are found from For exchange field in z-direction we have m 1 = m 2 = 0. The minimum can thus be found on the k z -axis by minimizing the same function f (c) Eq. (108) as for model I. We thus find for the critical value V cr,z in this case: From this expression we see that for model II in general the critical value V cr depends on the direction of the exchange field, in contrast to model I, where it was isotropic. The position of the Fermi points in this case is given by Eq. (112) for V > V cr,z . Again, one sees that the position of the Fermi points in model II varies with the direction of the exchange field, in contrast to model I.
To evaluate this expression analytically it is useful to employ the residue theorem in the following way [67]: for given k x and k z the quantity γ(k) = detD(k) = m 2 0 +m 2 1 +m 2 2 +m 2 3 −V 2 x −V 2 y +2im 1 V y −2im 2 V x (121) for k y ∈ [−π, π] defines a closed path in the complex plane. Therefore, by transforming z = detD(k) the winding number can be expressed as the integral w(k x , k z ) = 1 2π Im˛γ dz z = I(γ, 0).
where I(γ, 0) is the winding number of the curve γ around the origin. Thus, if the curve does not enclose the origin, the winding number becomes zero and no zero energy surface state exists. Now, the path γ crosses the real axis, when Im detD(k) = 0. This happens, when For a given k x this equation has two solutions k y,1 and k y,2 = π − k y,1 , if Vy Vx sin k x ∈ (−1, 1). The path encloses the origin only, if these two crossings possess different sign of the real part, i.e. the criterion for the existence of the flat band becomes 0 > γ(k x , k y,1 , k z ) γ(k x , k y,2 , k z ). (124) The boundary of the flat band in k x -k z -space is then given by the criterion 0 = γ(k x , k y,1 , k z ) γ(k x , k y,2 , k z ).
Now, as on the Fermi surface we have detD(k) = 0, it follows that the boundary of the flat band is just the projection of the (one-dimensional) Fermi surface onto the k x -k z -plane.