Impurity scattering in Weyl Semimetals and their stability classification

Weyl Semimetals (WS) are a new class of Dirac-type materials exhibiting a phase with bulk energy nodes and an associated vanishing density of states (DOS). We investigate the stability of this nodal DOS suppression in the presence of local impurities and consider whether or not such a suppression can be lifted by impurity-induced resonances. We find that while a scalar (chemical potential type) impurity can always induce a resonance at arbitrary energy and hence lift the DOS suppression at Dirac/Weyl nodes, for many other impurity types (e.g. magnetic or orbital-mixing), resonances are forbidden in a wide range of energy. We investigate a $4$-band tight-binding model of WS adapted from a physical heterostructure construction due to Burkov, Hook, and Balents, and represent a local impurity potential by a strength $g$ as well as a matrix structure $\Lambda$. A general framework is developed to analyze this resonance dichotomy and make connection with the phase shift picture in scattering theory, as well as to determine the relation between resonance energy and impurity strength $g$. A complete classification of impurities based on $\Lambda$, based on their effect on nodal DOS suppression, is tabulated. We also discuss the differences between continuum and lattice approaches.


Introduction
The history of relativistic (Dirac) fermions in solid state band structures has been known since Wallace [2], who first considered a single layer of hexagonal graphite, i.e. graphene. It was however generally believed that such structures were intrinsically unstable and impractical to fabricate, but advances in materials preparation and experimental techniques have led to a surge of interest [3] and a new class of materials, known as Dirac materials [4][5][6][7]. These materials have one or more symmetry-protected Dirac nodes where the density of state (DOS) vanishes and about which the energy dispersion is linear. The Dirac nodes are topological in nature, appearing as vortices or monopoles in the bulk Brillouin zone, and their presence typically requires some fine tuning such that the appearance of the Dirac structure marks a quantum phase transition between gapped phases. The properties of the gapped phases on either side of the transition differ in terms of their boundary behaviors, with the surface or edge spectra reflecting the topology of the bulk band structure. Familiar examples of systems with Dirac nodes include graphene and the class of materials known as topological insulators [3,5,6].
What distinguishes the Weyl semimetals (WS) within this framework is that the nodal structure exists not uniquely at a quantum critical point, but throughout an entire phase. In three-dimensional (3D), a Dirac point consists of two Weyl nodes of opposite chirality overlapping at the same point in k-space. While crystal symmetry may protect the two Weyl nodes from coupling [8,9], thereby stabilizing the Dirac node, under a general perturbation they will be coupled and thereby open up a gap. In a system invariant under both timereversal (T ) and inversion (I) symmetries, the nodal structure is at least four-fold degenerate (Kramer's degeneracy), and the aforementioned separation requires the breaking of either T or I symmetry, or both. This can be achieved for example by introducing external field, magnetic bulk impurities, electron interaction, etc [10][11][12][13]. If the perturbation which lifts the degeneracy is strong enough, a new gapless phase may result, which is the WS phase. The WS has been identified in several recent studies [17][18][19][20][21][22][23], and is a stable phase because further bulk perturbations can only shift the Weyl nodes without eliminating them. Recent work has also explored more exotic band structures involving Weyl fermions, for example their coexistence with quadratic (massive) fermions [14,15], or the symmetry-enforced overlap of Weyl nodes of the same chirality [16]. It is known that band insulators may undergo a topological phase transition as the bulk gap collapses and re-opens again; in this sense, the WS-whose bulk gaps are closed-is an intermediation of two topologically distinct gapped phases [24], e.g. from a trivial insulator to a topological insulator [17], or a Chern insulator if T is broken [10].
A minimal model of the WS was constructed by Burkov, Hook and Balents (BHB) [1]. BHB considered a massless 3 + 1 dimensional four-component Dirac fermion model, initially with time-reversal and inversion symmetries. A k · p expansion about the T and I -symmetric Brillouin zone center yields the four-component Hamiltonian H 0 (k) = 3 a=1 k a a + m 4 , written in terms of Dirac matrices. The sign of the mass m tells if the bulk insulating state is normal or topological. BHB showed that adding a homogeneous T and/or I breaking term initially gaps out the Dirac spectrum, but that for sufficiently large symmetry breaking a WS phase appears, with Weyl nodes occurring at two distinct k points. (In some cases, the two central bands touch along a circle in k-space. ) From the perspective of bulk-boundary correspondence, Weyl nodes may appear as the ends of the so-called Fermi arc [20], the locus of gapless surface states interpolating between the projections of Weyl nodes on the surface Brillouin zone. Such gapless modes participate in surface transport, with their multiplicity proportional to the arc length. This gives rise to an anomalous Hall effect of the T breaking WS, which recently has been shown to survive even when the Weyl nodes are subsequently gapped out by node-mixing scatterings, and is attributed to the persistence of chiral anomaly [25].
In this paper, we investigate the effects of localized impurities on the bulk electronic structure of the WS. In particular, we address the question of whether or not the DOS suppression at Weyl nodes in clean samples can be lifted via impurity scattering. Local impurities are modeled as V = g δ(x) where g ∈ R is the coupling strength, δ(x) restricts the impurity to the site at x = 0, and is a matrix structure encoding its physical type, e.g.
= I for scalar (chemical potential) impurity, and ∝ σ z for a magnetic impurity polarized along the z direction. We will speak of the stability of an energy ω under the scattering of a -type impurity in the following sense: if a resonance or bound state can be induced at ω by with some g, then ω is unstable with respect to . Otherwise it is stable. Close to the nodal energy, the DOS vanishes as ω 2 . If the nodal energy is unstable, the resulting resonances will give rise to sharp peaks in the DOS which disrupt the pristine Dirac spectrum [26]. Bulk transport consequences for scalar impurities were considered in [1,11,[27][28][29]. [28] studied the effect of rare regions in a dirty WS. Effects of scalar and magnetic impurities on the surface Dirac nodes of 3D topological insulators were studied in [30].
One might be tempted to draw intuition from the more familiar single-band problems and conclude that an impurity can induce resonance or bound state at arbitrary energy, given the freedom in choosing its strength g, making all energies unstable. We find that in the multipleband case such as the WS, while this still holds for scalar impurities, it is not true in general. Instead, stability depends crucially on the type of impurity, which is mathematically classified by its commutation relation with the matrices in the local Green's function. For some impurities, resonances and bound states are forbidden over a wide range of energies. Typically, an impurity is a foreign atom or local crystalline defect in an otherwise pristine material. Thus a realistic impurity potential should always involve a local scalar scattering component. If this scalar effect dominates the impurity, intragap resonances can be induced which will destabilize the Weyl node at a single particle level. If, on the other hand, the scattering is dominated by the resonance-forbidding components, then the Weyl node will remain. This paper is organized as follows: in section 2, we present a general framework to address the existence of impurity resonances, and the dependence of their energies on the impurity strength. In section 3, we introduce a four-band tight binding lattice model of the WS in terms of the matrices, adapted from the continuum BHB model [1]. Before analyzing the impurity effect in this lattice model, we first discuss in section 4 the situation in the low energy theory, namely the original continuum BHB model, and show that a natural momentum cutoff around the Weyl nodes in such theories dismisses the important physics of a stabilized Weyl node. We then turn to a full lattice treatment: in section 5, we apply the method developed in section 2 to the WS model and classify impurities according to their effect on the electronic structure. We shall throughout the paper restrict to being one of the sixteen matrices (including I); linear combinations thereof can be analyzed in the same fashion but with more tedious algebra. We will first illustrate the method in section 5.1 with the simpler case of a Dirac semimetal which is invariant under both T and I. The (fine-tuned) Dirac node is shown to be unstable with I even impurities, but stable with I odd ones. We apply the same approach to WS with a symmetry breaking term η where is a matrix and η is its strength. For a generic energy ω, we find that its stability depends on both the impurity type and the symmetry breaking term η . Section 5.2 discusses I breaking WS (which may or may not break T) where anticommutes with I. In this case, the Weyl node energy is found to be stable for any that does not commute with the local Green's function, but unstable if it commutes. Section 5.3 discusses the I -symmetric WS, in which necessarily breaks T. Again, impurities commuting with the local Green's function will disrupt the Weyl node stability. Those that do not fully commute yield either a nodal energy stable over the full range of η, or a critical symmetry breaking amplitude η c -which is fully determined by parameters of the clean system-that splits the η axis into two phases where the nodal energy is stable in one phase and unstable in the other. The critical amplitude η c is found to be related to a type of band inversion and indicates a phase transition in the gap of resonant impurity band structure, reminiscent of band inversions in topological/Chern insulators that are responsible for topological phase transitions. We conclude in section 6.

Resonance criteria in a generic multi-band system
The effect of localized impurities can be studied in general using the standard T -matrix formalism [31,32]. We briefly summarize the procedure below and establish notation. Given a Hamiltonian H = H 0 + V , its Green's function is where z ∈ C is the complex frequency, G 0 (z) ≡ (z − H 0 ) −1 is the Green's function of H 0 and T = V (I − G 0 V ) −1 is the T -matrix. Assume the impurity potential V is localized at a spatial point r = 0: V r r = g δ r,0 δ r ,0 , where g is the potential strength and is a dimensionless matrix whose rank is equal to the number of bands. The Green's function connecting r and r is G r r = G 0 r r + G 0 r0 T 00 G 0 0r , and T 00 (z) = [g −1 −1 − G 0 00 (z)] −1 is the only nontrivial block of the T -matrix. For translationally invariant systems, the local Green's function is where H 0 k is the Fourier transform of H 0 and N is the number of k points, i.e. the number of unit cells of the crystal. The corresponding local DOS (LDOS) at site r with energy ω is ρ r (ω) = −ImTr G r r (ω + i 0 + ).
Bound states and resonances are consequences of the energy spectrum reconstruction induced by impurities. Before taking the thermodynamic limit, eigenvalues of H are poles of T (ω) on the real ω-axis, i.e. the zeros of T −1 (ω) = g −1 −1 − G 0 00 (ω). Upon tuning of the impurity strength g −1 , each pole will trace out a curve in the (ω, g −1 ) plane. Note that g −1 = ±∞ are identified as g = 0 and ±∞ as g −1 = 0. Thus one may consider an adiabatic cycle in which g −1 : −∞ → 0 → +∞. After a full cycle, the poles must collectively recover their initial positions, which are the set of clean states at g −1 = ±∞. Individual poles may either stick close to one clean state, or migrate between different ones. An example is shown in figures 1(a) and (b) for the spectral evolution of a graphene sheet of 6 × 6 and 10 × 10 unit cells, respectively, in the presence of a scalar impurity.
A point (ω, g −1 ) on the spectral evolution curves outside the bands of the clean system represents a bound state of energy ω induced by an impurity of strength g. For those inside the clean bands, one has to distinguish between poles very close to-hence mere perturbations of-a clean state, and those in the middle of a migration. The latter, like bound states, are manifestations of the impurity effect. They differ in that an increase in system size N has little influence on the bound states, but will split spectral lines inside the clean bands to accommodate newly created clean states; see figures 1(a) and (b) for this lattice size effect. What remains unchanged when increasing N is the trend of rapid pole migration.
An effective way to extract the locus of such rapid migrations, which we identify as resonances, is to find the zeros of is the Hermitian part of the retarded local Green's function, where the imaginary part is taken to be greater than the energy spacing between consecutive bulk levels. In so doing, divergences in T originating from poles of G 0 00 (ω) itself are eliminated from the zeros of T −1 (ω), leaving only those caused by the aforementioned spectral migrations. The zeros of T −1 (ω) are shown for the case of graphene as green curves in figures 1(a) and (b). In these figures, they sit close to the inflection points of the spectral curves, where one might say the pole migration is most rapid.
In the thermodynamic limit N → ∞, the Hermitian matrix T −1 remains well defined. Since → 0 + , T −1 and T −1 are identical outside bulk bands, the zeros of T −1 can be used to identify both bound states and resonances. Callaway [33] showed that the phase shift at ω, viz δ(ω) ≡ arg det T (ω + i 0 + ), equals π × N (ω) where N (ω) is the difference in the total number of states of H and of H 0 below ω. One definition of resonance in this context is for the phase shift to be ±π/2, viz. Re det T = 0, the reason being that the number of extra states is a half-odd-integer, which represents the 'center' of the process in which one extra state is gained or lost. The former is called a resonance and the latter an anti-resonance. Over the full range of ω they must balance each other if introduction of an impurity does not change the total number of states, a form of the Friedel sum rule. This is related to our resonance criteria in that the zero of T −1 is the center of a spectral migration-a spectral line migrating past ω, by definition, contributes | N (ω)| = 1 to the total number of states below ω. The difference between our criteria and the phase shift picture is that there are in general N B branches of spectral evolutions at any ω where N B is the number of bands, e.g. N B = 2 in the graphene example shown in figure 1. where γ k = 1 + e −ik 1 + e −ik 2 and k i = k · a i , with a i (i = 1, 2) being the two primitive direct lattice vectors. Panels (a) and (b) show the spectral evolution of H = H 0 + gI (solid and dotted gray curves) with impurity strength g, for different lattice sizes. These are solved as zeros of the eigenvalues of the , with the eigenvalue index encoded by different line color/types. Migration of such levels between different clean states (g −1 = ±∞ limit) constitutes a resonance (inside bulk bands) or a bound state (outside bulk bands), and can be extracted as zeros of T −1 (ω, g), the Hermitian part of T −1 , shown as dashed green curves, with solid circular blue points overlaying on the branch corresponding to the solid spectral flow and red empty square on the branch corresponding to the dotted spectral flow (see the text). The discontinuity in the green curves at ω = ±1 is concomitant with the Van Hove singularity in the DOS (not plotted) (c) shows the phase shift arg det T (ω + i 0 + , g) in the thermodynamic limit, where ± π 2 (heaviest red/blue) could be interpreted as resonance or bound state, see the text. (d) Plots the norm of (Continued) retarded T matrix, which can be used to distinguish between resonance and anti-resonance that is hard to tell from (c), the former corresponding to the dark feature and the latter suppressed in such a plot. The dashed green curves in (c) and (d) are the same as those in (a) and (b).
Our criteria is essentially to track the increment/decrement contributed by any single branch, whereas Re det T = 0 takes into account all branches. In any case, the difference is consequential only in identifying the location of the anti-resonances. For comparison, we plot the phase shift (color map) together with the zeros of T −1 (dotted green line) for a graphene impurity in figure 1(c).
The distinction between resonance and anti-resonance is not immediately apparent from the phase shift plot. It relies on the sign of s = ∂δ/∂ω: s > 0 is a resonance and s < 0 an antiresonance. Furthermore, in cases where two bound states/resonances are close together-or even degenerate as can happen in Dirac semimetals to be discussed in section 5.1-the phase shift will experience a 2π change over a small ω window, which is equivalent to zero numerically and hence hard to resolve. A more transparent way is to plot the matrix norm ||T 00 (ω || is large, then the DOS is in general enhanced, and one obtains a resonance. This is shown for the graphene example in figure 1 To facilitate analytical treatment, we will henceforth use the T −1 approach and make no distinctions among resonance, anti-resonance and bound state. All of them will simply be referred to as 'resonance'. As needed we will also exhibit ||T (ω + i 0 + )|| plots, where antiresonances are suppressed, as a check.
Consider next the existence of a resonance at an arbitrary point in (ω, g −1 ) space, i.e. the condition for at least one eigenvalue of T −1 (ω, g −1 ) to be zero. If g is allowed to be complex, then there are as many solutions of g at a given ω as the number of bands: g −1 = u a (ω), where {u a (ω)} are the eigenvalues of G 0 00 (ω) . However, only real g is physical. Thus existence of resonance demands at least one eigenvalue of G 0 00 to be real. The required coupling strength is g = 1/u a (ω). An immediate corollary is that impurities with [ , G 0 00 ] = 0 can induce a resonance at arbitrary energy for some g, because product of commuting Hermitian matrices has real eigenvalues. Single-band problems fall in this category ( = I).

Weyl semimetal (WS) models
We now apply the method described above to lattice systems adapted from the continuum models of BHB [1]. The k-space Hamiltonian in the -matrix basis is where the k-dependent coefficients are taken to be and η is k-independent for simplicity. is a matrix which breaks time-reversal (T ) and/or inversion (I ) symmetry to be defined later. The following matrix convention is used: where τ i and σ i are two sets of Pauli matrices acting on the orbital and spin degrees of freedom, respectively. In this model, t and t are intra-orbital hoppings, while t 1 is the (spin-mixing) hopping between different orbitals. Different conventions of matrices may have different physical interpretations: for example, in the convention above, 12 = −I ⊗ σ z represent a magnetic field in the z direction, but in other conventions it may also have orbital effects. Results obtained below are independent of the convention used.
BHB showed how the emergence of stable point or line nodes in the spectrum of H 0 beyond a critical perturbation strength depends on the symmetry of the perturbation [1]. Define timereversal as T = KR where K is complex conjugation and R = I ⊗ i σ y , and inversion as I = 4 . Symmetry properties of all matrices can be found in table 1. If η = 0, then the model is both T and I -symmetric. In this case, fine-tuning λ = 0 creates a Dirac node at k = 0 where the four bands converge to the energy E = −6t − ε 0 . While the Dirac point is gapped out by nonzero λ, and is therefore unstable in the T and I -symmetric system, it splits into two Weyl nodes if either T or I is broken by [1,17]. In this case the nodal structure survives in a range of parameters and constitutes a stable nodal phase, i.e. the WS phase.
In BHB's language, H 0 with η = 0 is the unperturbed Hamiltonian, and the η term is the symmetry-breaking perturbation. Although the Weyl nodes and hence the semimetal phase are stable under such homogeneous bulk perturbations, the characteristic suppression of DOS at the nodal energy in the WS phase, may be destroyed by another type of perturbation-namely localized impurity potentials-if resonances can be induced at the nodal energy. Hereafter, we shall refer to the full H 0 of equation (2), including the homogeneous term η , as the unperturbed Hamiltonian, and regard the local impurity potential as the perturbation.
The model of equation (2) is lattice-based, which entails a specific cutoff structure, and hence not generic like the BHB Hamiltonian H BHB = 3 i=1 k i i + m 4 + η . It will be instructive to first look at the effect of scattering from the more universal low energy states living in the vicinity of Weyl nodes, which we shall analyze in section 4 for the BHB model. However, as we shall see there, the very act of taking a momentum cutoff will leave out the possibility of a stable Weyl node. This is not surprising because spatially localized impurities are homogeneous in the momentum space and inevitably scatter high momentum states. Thus a lattice treatment is necessary. For the lattice theory, we will first discuss the impurity effect in the T and I -symmetric system (η = 0). While it does not yield a WS phase, it is simple enough to be used as a demonstration of the general framework outlined in section 2. Then we will move on to unperturbed systems with T and/or I broken (η = 0) where a WS phase does exist. Since the Hermitian part of the local unperturbed Green's function, G 0 00 , is of central importance to the impurity classification, we will first classify the unperturbed system according to the type of matrices appearing in G 0 00 , and then for each of them, classify the impurities according to their commutation with G 0 00 .

Impurity scattering in the low energy theory
In this section we focus on the low energy theory described by the following BHB Hamiltonian: In both cases the scale of nodal momentum is = η 2 − m 2 . For η > m > 0, = 21 generates two point nodes at k = (0, 0, ± ), whereas = 35 generates a nodal line at k = ( cos φ, sin φ, 0). In the vicinity of the Weyl nodes, the s = −1 bands have linear dispersion, whereas the s = +1 pair is gapped, If we are looking for low energy resonances, ω ∼ 0, then since the Green's function is weighted by 1/(ω − E), it is reasonable to focus on the scattering of low energy states by (i) projecting onto the s = −1 bands and (ii) adopting a momentum cutoff such that only momenta within a distance Q from the nodes are considered-a sphere around point nodes and a tube around nodal line-with ω Q η. In this approximation the local Green's function for both s have the following form: see appendix B for derivation and equation (B.54) for the expression of R(m/η). The Hermitian G 0 00 (ω) is obtained by taking the real part of a(ω).
Let us now analyze the resonance condition. The T −1 matrix is where a r is the real part of a. The impurity either commutes or anti-commutes with 4 since the latter is itself one of the 16 matrices. If they commute, then det T −1 00 = 0 yields real solutions for g −1 , i.e. resonance could be induced and ω ∼ 0 is unstable. If, on the other hand, and 4 anti-commute, then det T −1 00 = 0 gives Now, by m 2 + 2 = η 2 , one obtains g −1 = ±|a r (ω)| /η, which is still real, i.e. ω ∼ 0 is unstable for anti-commuting . Thus in the cutoff scheme adopted here, ω ∼ 0 is unstable regardless of the type of impurity. Note however that in the case of anti-commuting ones, the impurity strength g always comes in ± pairs, whereas in the commuting case it might not.
An essential difference between the commuting and anti-commuting impurities is that in the latter case, the solution for g −1 contains a square root. When higher momentum states are considered, the argument of the square root may become negative, and stabilize the nodal energy.
To see this, note that without the low energy restriction, the local Green's function of the BHB Hamiltonian with = 21 and 35 has the form (see equation (A.20) where the coefficients are in which · · · denotes the k-space integral (for continuum) or sum (for lattice), and In the low energy approximation, b 1 and b 2 vanish due to the momentum cutoff. This is because at the Weyl nodes, E 2 + = d 2 ⊥ = 4η 2 , thus in the low energy approximation, the second and third terms inside · · · of both b 1 and b 2 cancel, yielding where V cutoff /V BZ is the ratio between the volume within the momentum cutoff and that of the first Brillouin zone. This volume ratio comes from the evaluation of 1 , and is of the order (Q/2π ) codim ∼ (Q/η) codim → 0 in the low energy approximation, with codim being the codimension of the Weyl node, which is 3 for a point node and 2 for a nodal line. When higher momentum states are included, b 1 and b 2 will no longer be suppressed, and will change the argument under the square root, possibly making it negative and stabilizing the nodal energy. A more careful analysis necessitates a lattice treatment, which is what we shall do in the rest of the paper.

T and I -symmetric H 0 (k)
We now return to the lattice model of section 3. To illustrate the resonance criteria of section 2, we first consider the case with η = 0. Inverting equation (2) yields the unperturbed k-space Green's function The i (i = 1, 2, 3) terms will vanish after summation over k due to the oddness of d i (k), as required by inversion symmetry. The local Green's function is thus where N is the number of k points, and As discussed in section 2, the existence of a resonance depends on whether or not the eigenvalues of T −1 00 , i.e. the Hermitian part of the inverse local T -matrix, can be zero. Since the only matrix in the G 0 00 decomposition is 4 = I, there are only two classes of impurities according to their inversion property.

Inversion-even impurity.
In this class we have = I, 4 or µν (µ, ν = 1, 2, 3, 5). Since they all commute with G 0 00 , a resonance can be induced at arbitrary energy, i.e. a solution of det T −1 00 = 0 exists for real g. To illustrate, we solve for g(ω), the value of g which produces a resonance at energy ω, for all three cases:

1.
= I : this is a scalar impurity, and The principal values of a and b are implicitly taken. Setting the lhs to zero yields These are shown as the light blue dashed lines in figure 2(b).

2.
= 4 : this impurity flips the sign of the inversion-odd component, yielding The resonance condition is thus
(28) The results are plotted as dashed lines in figure 2(a). Thus inversion-odd impurities cannot induce any resonance in the range of energy ω for which |a(ω)| < |b(ω)|.

Band center approximation.
To understand the general trend of the resonance solutions g(ω) and get a sense of the stability region, it is useful to obtain an approximation for the expansion coefficients a(ω) and b(ω). To this end we introduce the band center approximation (BCA): a generic k-space Hamiltonian H (k) can be written as H (k) = H 00 + δ H k where H 00 = H (k) is the local Hamiltonian and · · · denotes k-space averaging. The eigenvalues of H 00 can be thought of as some sort of average energy of the bands of H (k) (band centers). LetḠ then the local Green's function is where we have used δ H k = 0. The BCA amounts to replacing the local Green's function, G 00 , with the Green's function of the local Hamiltonian,Ḡ. Equation (31) is an expansion in powers of δ H k /(ω − H 00 ), where the numerator is roughly the band width, and the denominator is the distance from ω to the band centers. The BCA works well if the distance of ω from some band center, say that of band A, is greater than A's bandwidth. Note that such an ω, although outside band A, may well be inside another band, say band B. From the BCA point of view, a resonance at some ω inside band B is actually a consequence of the coherent superposition of states mainly in some other band (A). The multiple-band scenario is to the benefit of the BCA.
Applying BCA to the η = 0 model here, one finds from equation (2) that H 00 = αI + β 4 where α = −ε 0 and β = −12t − λ, hencē G and G 00 have the same form of decomposition. This will prove useful in the more complicated situations where is present. For the scalar potential = I, the BCA resonance solution g −1 ā(ω) ±b(ω) = [ω − (α ± β)] −1 resembles two hyperbolae centered around the band centers α ± β. They can be identified qualitatively from the light blue dashed curves in figure 2(b), although numerically the two branches of each hyperbola, instead of being divergent, are connected around their respective band centers due to higher order effects in equation (31). For I odd impurities, such as the magnetic impurity = 12 , the stability condition equation (29) implies α − |β| < ω < α + |β|, i.e. stable energy ω is bounded by the two band centers, as can be seen from figure 2(a). This region in particular includes the (fine-tuned) Dirac point or the central gap. We thus conclude that the Dirac node is generically stable for I odd impurities and unstable for I even impurities.

H 0 (k) with I breaking
We now consider the case where breaks inversion. In the BHB scheme of alternating layers of 2D topological and nontopological insulators (T1/N1) stacked along say the z direction, this can be realized for example by applying a voltage bias across each TI layer, breaking the inversion symmetry between the top and bottom surfaces of each T1 layer. According to table 1, = µ4 (T even) or µ (T odd) where µ = 1, 2, 3, 5. The local Green's function has the decomposition While this decomposition can be obtained analytically (see equations (A.14), (A.15), (A.17) and (A.18)), its structure is easier to understand from the BCA: symmetry consideration demands that the local Hamiltonian H 0 00 ≡ 1 N k H 0 (k) = α I + β 4 + η where the first two terms are the only possibilities to conserve both T and I, see table 1. Its inverse can potentially have four terms, I, 4 , and 4 . Since anti-commutes with 4 , their cross term must vanish, yielding the form in equation (33). We will come back to BCA later.
The resonance condition can now be solved for different impurities. As an example, consider = µ4 and = µμ = −1 whereμ ∈ {1, 2, 3, 5} \ {µ}. This type commutes with 4 but anti-commutes with , and includes the purely magnetic impurities 12 , 23 and 13 . The Hermitian part of the inverse T 00 matrix is Using the anti-commutation {g −1 µμ − b 1 4 , µ4 } = 0, the above can be rearranged into Setting T −1 00 = 0, both sides can be simultaneously diagonalized, and the eigenvalues of the rhs are (g −1 ± b 1 ) 2 . The condition for vanishing T −1 00 is thus Resonance then requires g −1 to be real, viz. |a(ω)| > |b 2 (ω)|. The occurrence of a possibly negative term under the square root stems from the anti-commutation of with , i.e. the interplay between the bulk symmetry breaking field and the impurity. Similar analysis can be carried out when is any of the sixteen matrices. The results are summarized in the third column of table 2. The sixteen -matrix impurity candidates can be classified into four classes labeled by their commutation with 4 and µ4 : (+, −) denotes commuting with 4 and anti-commuting with 4µ , and similarly for (+, +), (−, +) and (−, −). Impurities belonging to the same class have the same resonance condition. A nontrivial solution arises if there is at least one anti-commutation, giving rise to a possibly negative term under the square root, and the protection of DOS suppression at Weyl nodes. The unperturbed H 0 is parameterized by the symmetry-breaking strength η, and one can ask how it affects the system's ability to induce resonance at energy ω. The (ω, η) space is thus divided into two phases according to the existence of resonance. These are shown in figure 3, where the stable phases (no resonance) are colored.
The shape of the phase boundaries can be qualitatively understood in terms of the parameters of the Hamiltonian using BCA: the local Hamiltonian is H 0 Note that these coefficients are independent of , thus the stability of the Weyl nodes can be predicted according to the impurity class, regardless of . The BCA version of the phase boundaries are shown as dotted lines in figure 3. The DOS suppression at the bulk Weyl nodes is protected for impurities in classes (+, −), (−, +) and (−, −). The only unstable class is (+, +), due to its fully-commuting nature with G 0 00 . Table 2. Impurity classification for I breaking Weyl material. The symmetry breaking term in H 0 (k) is = µ4 (T even) or µ (T odd) where µ ∈ {1, 2, 3, 5}. The second column indicates commutation (+) or anti-commutation (−) of the impurity matrix with 4 and , respectively. Elements in each class are enumerated in the first column: if the cell has two sub-cells, the left one corresponds to = µ4 and right one = µ ; otherwise the enumeration is identical for both . The value of g at which T −1 00 has a zero eigenvalue is listed in the third column, and the condition for it to be real (the resonance condition) is shown in the fourth column. The fifth column shows the resonance condition as given by the BCA, which are simple expressions in terms of the Hamiltonian parameters. Note that the values of the Green's function coefficients a(ω) and b i (ω) depend on the choice of that breaks I, but the BCA conditions are independent of . The stability of Weyl nodes, if they exist, is listed in the last column.
Any ω Any ω Unstable

H 0 (k) with I -symmetric
To split the Dirac node into two Weyl nodes, an I -symmetric must break T. Thus = µν (µ = ν = 4) according to table 1. The local Green's function is Note that all three matrices in equation (36) mutually commute, and the product of any two is equal to the third. This implies that the impurity either commutes with all of them, or it commutes with one and anti-commutes with the other two (because the product of any two commutation signs should produce the third). In the fully commuting case, resonance can where s 4 , s µν and s are eigenvalues of 4 , µν and respectively and independently take the values ±1. For other , there are two anti-commutations. We can relabel the three matrices in equation (36) according to their commutation with , and write the inverse T matrix as where The two parentheses in equation (38) mutually commute, thus T −1 00 is block-diagonal: in the eigen-subspace of C with eigenvalue ±1, the matrices C , , A and A reduce to 2 × 2 blocks, denoted as ±I, ± , ± A and ± A , respectively, all of which square to I. Since the projectors onto the subspaces of C commute with C , A , A and , the mutual (anti)commutation relations of the latter four are inherited in both subspaces. Setting T −1 00 = 0 in equation (38) for both blocks then yields Squaring both sides and using the fact that ± A ± A = ±I, which follows from A A = C , one gets The resonance condition is for g to be real, viz. Table 3. Impurity classification for I -symmetric (hence T breaking) Weyl material. The class to which belongs are labeled by the three signs of the commutation of with 4 , µν , and 4 µν in that order, where + denotes commute and − anti-commute. The two indices µ, ν ∈ {1, 2, 3, 5} and are fixed by the unperturbed Hamiltonian. The index p ∈ {1, 2, 3, 5}\{µ, ν}. s and s take the values of ±1. Solution of g yielding det T −1 = 0 are summarized in the third column, see equations (37) and (41) in text. The resonance conditions for each class can be deduced by requiring g −2 to be positive (so that g is real), and are explicitly spelled out in the fourth column using BCA, which only need to be satisfied for either s = 1 or −1. Stability of the Weyl nodes are listed in the last column. α = −ε 0 and β = −12t − λ.
Class g −2 Resonance (BCA) Stability  if at least one of ± is satisfied. This is enumerated in the third column in table 3 and plotted in figure 4. As before, the expansion coefficients in equation (36) can be estimated by the BCA and used to approximate the boundaries between the resonant and nonresonant phases in the (ω, η) space. The local Hamiltonian is H 0 00 = αI + β 4 + η µν where α = −ε 0 , β = −12t − λ. Its Green's function isḠ(ω) =ā(ω)I +b 1 (ω) 4 +b 2 (ω) µν +b 3 (ω) 4 µν with where The resulting resonance conditions are summarized in the fourth column of table 3, and the conditions for stable Weyl nodes (if exist) in the last column.
The stable zones of the two classes (+, −, −) (green in figure 4) and (−, +, −) (red in figure 4) are restricted to opposite sides of a critical value of the symmetry-breaking strength η = η c . Furthermore, near η c , the region of stable energy narrows down toward the Weyl node. One can think of the resonance energy as forming an impurity band generated by a continuum of impurity strengths g. Then the stable zones constitute gaps in such bands. In this sense, η c marks a phase transition of the impurity band from gapless to gapped. The existence of η c can be understood from the BCA, according to which the phase boundaries are given by shown as gray dotted lines in figure 4. These are the two central band centers (eigenvalues of H 0 00 ). They cross at η = |β|, which gives the critical strength η c . η c = 2 in figure 4. This is reminiscent of the bulk band inversion in topological/Chern insulators that signifies a gapless to gapped transition in their surface/edge spectrum.
To illustrate the above impurity band phase transition, we employ the average T -matrix approximation (ATA) to investigate the effect of an ensemble of local impurities spatially uniformly distributed with concentration c [31,32]. The entire ensemble has the same matrix form, but with strength g given by some distribution f (g). In the ATA formalism, statistical averaging over the f (g) will restore translational symmetry. The impurity effect is then captured by a local self-energy, loc = c T 00 [1 + cG 0 00 T 00 ] −1 , where · · · denotes the f (g) averaging. The self-energy-corrected local Green's function is G loc (z) = 1 N k 1/(z − H 0 (k) − loc ) and the average LDOS is ρ(ω) = −ImTr G loc (ω + i 0 + ). One then expect the ATA DOS, ρ ATA , to be enhanced from the clean fraction, (1 − c)ρ clean , for ω in the unstable phase but reduced in the stable phase (since the integrated DOS is conserved). This is shown in figure 5, in which we plot at a fixed η the simplest case where f (g) is a constant for g ∈ (0, 10] and zero otherwise. For this particular η value, the Weyl nodal energy is stable in (b) and (d), but is unstable in (a) and (c).

Summary and discussion
In this paper, we study the effect of localized impurities V = g δ(x) on the bulk electronic structure of WS. A general method is devised to detect whether or not a resonance can be induced at energy ω by a -type impurity. If such a resonance is possible, then ω is said to be unstable with respect to , otherwise it is stable. The stability of ω requires all eigenvalues of G 0 00 (ω) to have finite imaginary part. Here, G 0 00 (ω) is the Hermitian part of the retarded local Green's function G 0 00 (ω + i 0 + ). Otherwise, one can always use a coupling strength g = 1/u a (ω) to induce a resonance at ω, with u a (ω) being the purely real eigenvalue of G 0 00 (ω) indexed by a. The existence of real u a (ω) is equivalent to requiring the Hermitian part of the inverse T matrix to have a zero eigenvalue. An immediate corollary is that impurities commuting with G 0 00 can induce a resonance at an arbitrary energy, simply by tuning the impurity strength. This includes the physically important class of local chemical potential perturbations.
We applied this method to four-band lattice WS models, expressed in terms of Dirac matrices. For these models, the T matrix and its eigenvalues can be obtained analytically. Mathematically, one first classifies the clean WS according to whether or not inversion (I ) is broken. The difference is in the decomposition of their local Green's functions G 0 00 . In each case, impurities are then classified by their commutations with the matrices appearing in G 0 00 : anti-commutation with components of G 0 00 result in a square-root structure, which constrains the reality of u a (ω). Note that in this scheme, it is more relevant to know which matrices appear than their exact numerical coefficients. For this purpose the BCA-which replaces G 0 00 (ω) by (ω − H 0 00 ) −1 where H 0 00 is the local Hamiltonian-is quite useful as it has the same form of matrix decomposition with coefficients whose meanings are physically more transparent. Results for I breaking WS are reported in table 2, and for I invariant WS in table 3.
Realistic impurities are more likely to be linear combinations of multiple matrices, mixing orbital, magnetic, and chemical potential effects. A linear combination of impurities in the same class, if it still squares to identity, is no different from a single matrix in that class, and results obtained before hold unchanged. While other combinations are not studied here, it is reasonable to expect that stability will resemble that of the dominant component if there is one, and crossover will happen as the relative strengths change. We have confirmed this for several tractable cases of Dirac semimetals. The method for obtaining the relation between impurity strength g and the induced resonance/bound state energy ω may also prove useful in device engineering where specific energy levels are desired. For Dirac materials with random strength disorder, results similar to those shown in figure 5 are expected, where roughly speaking impurity induced states form their own band superimposed on the clean DOS, and stable energies constitute the band gap. Such impurity bands may modify transport properties if certain impurity 'superlattice' is approximately formed, or if the coherent length of singleimpurity resonances become compatible with the impurity density.