Rare region induced avoided quantum criticality in disordered three-dimensional Dirac and Weyl semimetals

We numerically study the effect of short ranged potential disorder on massless noninteracting three-dimensional Dirac and Weyl fermions, with a focus on the question of the proposed quantum critical point separating the semimetal and diffusive metal phases. We determine the properties of the eigenstates of the disordered Dirac Hamiltonian ($H$) and exactly calculate the density of states (DOS) near zero energy, using a combination of Lanczos on $H^2$ and the kernel polynomial method on $H$. We establish the existence of two distinct types of low energy eigenstates contributing to the disordered density of states in the weak disorder semimetal regime. These are (i) typical eigenstates that are well described by linearly dispersing perturbatively dressed Dirac states, and (ii) nonperturbative rare eigenstates that are weakly-dispersive and quasi-localized in the real space regions with the largest (and rarest) local random potential. Using twisted boundary conditions, we are able to systematically find and study these two types of eigenstates. We find that the Dirac states contribute low energy peaks in the finite-size DOS that arise from the clean eigenstates which shift and broaden in the presence of disorder. On the other hand, we establish that the rare quasi-localized eigenstates contribute a nonzero background DOS which is only weakly energy-dependent near zero energy and is exponentially small at weak disorder. We find that the expected semimetal to diffusive metal quantum critical point is converted to an {\it avoided} quantum criticality that is"rounded out"by nonperturbative effects, with no signs of any singular behavior in the DOS near the Dirac energy. We discuss the implications of our results for disordered Dirac and Weyl semimetals, and reconcile the large body of existing numerical work showing quantum criticality with the existence of the rare region effects.


I. INTRODUCTION
Recently, there has been an intense experimental effort to find gapless semiconductors that host isolated points in momentum space with linearly touching valence and conduction bands. This thrust has been fueled by the exciting possibility of studying massless three-dimensional Dirac (for Kramers degenerate bands) and Weyl (for non-Kramers degenerate bands) fermions in solid state systems. (The fact that the two-dimensional version of a Dirac-Weyl system already exists in the form of graphene has obviously been a great impetus in this search for three-dimensional Dirac-Weyl materials.) This has led to establishing three-dimensional Dirac semimetals in the compounds Cd 3 As 2 (Refs. 1-3), Na 3 Bi (Refs. 4 and 5), Bi 1−x Sb x (Refs. [6][7][8], BiTl(S 1−δ Se δ ) 2 (Refs. 9 and 10), (Bi 1−x In x ) 2 Se 3 (Refs. 11 and 12), and Pb 1−x Sn x Te (Refs. [13][14][15]. While, even more recently the existence of Weyl semimetals [16][17][18] in TaAs (Refs. 19 and 20) and NbAs (Ref. 21) has been established. This low energy description is also applicable to various other physical systems that host gapless Dirac or Weyl points such as the pyrochlore iridates 16 and the Bugliobov quasiparticle properties of nodal superconductors. With the experimental discovery of such a large number of Dirac-Weyl materials (and the great deal of interest and excitement surrounding them), as established by their electronic band structures through photoemission spectroscopy (i.e. linearly touching conduction and valence bands), one of the immediate important questions is how robust this noninteracting clean system is to the presence of interaction and disorder, physical effects invariably present in real solid state materials. Here we study the fundamental effects of static potential disorder on noninteracting Dirac-Weyl systems. (We note that typically these materials are considered to be weakly interacting due to the strong screening provided by the large background lattice dielectric constant in the systems.) Due to the invariable presence of disorder in all solid state materials there has been a substantial amount of theoretical activity studying the effect of disorder on noninteracting Dirac and Weyl fermions [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39] . Focusing on the undoped (i.e., Fermi energy at E = 0) Dirac point (i.e. the band touching point), the quadratically vanishing density of states at zero energy (ρ(E) ∼ E 2 ) associated with the linear three-dimensional energy band dispersion places these problems in a different class than that of a conventional metal with a parabolic energy dispersion and a nonzero Fermi energy. In a standard metal, the nonzero density of states at the Fermi level gives a finite mean free path at leading order in a random potential. (We note that a regular metal is different from a Dirac-Weyl system even in the hypothetical limit of a vanishing Fermi energy since there is an energy gap between conduction and valence bands for the regular metal whereas a Dirac-Weyl system is gapless-i.e. a regular metal simply becomes an ordinary gapped semiconductor in the zero Fermi energy limit whereas the gapless Dirac-Weyl system is a nontrivial semimetal for zero Fermi energy.) For the Dirac problem of interest here, from a scaling analysis of the action it is straightforward to see the perturbative irrelevance in three dimensions of disorder for massless Dirac and Weyl fermions 22 . Thus the semimetal (SM) phase could be stable up to some non-zero critical disorder strength, with a disorderdriven itinerant quantum critical point (QCP) into a socalled diffusive metal (DM) phase at higher disorder. (We mention here for completeness that in two-dimensions, e.g. graphene, disorder is perturbatively relevant, and thus infinitesimal disorder immediately converts undoped graphene from being a semimetal in the clean limit into a diffusive metal -thus two is the perturbative lower critical dimensionality for the disordered Dirac-Weyl problem.) The natural question that now arises with respect to the disordered three-dimensional Dirac-Weyl systems is whether the perturbative robustness of the semimetallic phase to disorder applies generally or is simply a perturbative result (perhaps to all orders in the perturbation theory), not valid in the nonperturbative theory. The goal of the current work is to settle this question definitively. Although the disordered Dirac-Weyl systems have been theoretically studied very extensively in the literature [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39] , essentially all of this work, except for a very recent one (Ref. 27), study the properties of the disorderdriven SM-DM quantum phase transition, taking it for granted that such a disorder-induced QCP indeed exists in three dimensions following the predictions of the perturbative field theory. Our current work reconciles the huge body of QCP theoretical work in the literature with the existence of nonperturbative or rare-region effects which lead to the 'suppression' or 'avoidance' of such a SM-DM QCP.
It is known that non-perturbative effects of rare regions may give rise to a non-zero (albeit exponentially small) density of states at zero energy for an infinitesimal strength of disorder 27,40 , thus converting the ballistic excitations in a weakly disordered SM to diffusive in the low energy limit, which thus results in the absence of a strict SM phase (with vanishing zero energy DOS). It is not uncommon for disorder to fundamentally change the nature of clean critical points (e.g. the Harris criterion 41 says this happens when the clean correlation length exponent ν < 2/d ), while the Chayes-Chayes-Fisher-Spencer (CCFS) 42 inequality for the exact correlation length exponent of the disordered system (ν ≥ 2/d) applies to critical points that occur in the presence of quenched randomness. Interestingly, the one loop perturbative renormalization group (RG) calculation of the critical exponents for the proposed SM to DM QCP are consistent with the CCFS inequality (since ν = 1, Refs. 22 and 23) as, in fact, are the 2-loop RG calculations 26,39 and all numerical estimates in the literature 24,25,32,33,35,36 , and therefore it is not a priori obvious that rare region effects should change the universality of this transition. Given the field theoretic RG analyses and the large body of direct numerical studies of the disorder-driven SM-DM QCP finding the various critical exponents and identifying the critical coupling as well as the apparent consistency between the theoretical (and numerical) correlation exponent with the CCFS inequality, it seems reasonable to assume that the rare regions arising out of nonperturbative disorder effects do not change the nature of the QCP in any substantial manner. In this work we explore the fate of this SM-DM QCP using specialized numerical techniques that allow for the direct study of these rare region effects. Our work definitively establishes that the putative SM-DM QCP becomes avoided or hidden due to the rare region effects, although the crossover effects of the avoided QCP show up in the numerical results. This is thus consistent with all the earlier numerical work finding an apparent existence of the QCP, which, we now argue, strictly speaking does not exist when examined to the lowest energies. Our work leads to the important conclusion that there is no disorder-driven SM-DM QCP in three-dimensional Dirac-Weyl systems, only an avoided QCP, and the Dirac point develops nonzero (albeit exponentially small) DOS even for weak disorder.
To put the problem into context, we first review the existing evidence for the disorder-driven SM to DM QCP in undoped Dirac-Weyl systems. The seminal work of Fradkin 22 established the existence of this disorder driven (perturbatively accessible) QCP. More recently, the properties of this proposed QCP have been calculated in a renormalization group treatment of the problem 23 , which has now been extended to two loops 26,39 . The field theory of the QCP can be constructed in terms of an interacting "Q 4 theory" (similar to a φ 4 theory for magnetism but now Q is the replicated matrix field) strongly coupled to massless Dirac fermions 33 , while for the Weyl case due to topological considerations a separate field theory has been derived 29,38 . Tuning the clean model away from the Dirac/Weyl limit by varying the power law of the dispersion relation 30,31 , this transition has been shown to occur even in some one-dimensional models 34 (akin to a longrange Ising model). Thus, the existence of the putative SM-DM QCP seems to be well-established from a field theoretic perspective.
Due to the non-interacting nature of the problem, various numerical techniques which can reach rather large system sizes, have been used to study the properties of this QCP. A main focus has been the direct calculation of the low energy density of states (DOS) ρ(E), since the DOS is expected to be singular at E = 0 at the transition. Moreover, the DOS can be related to the critical exponents via the scaling hypothesis 24 . Following this, the dynamic exponent z and the correlation length exponent ν have now been numerically estimated for several models using the directly numerically calculated DOS 24,[33][34][35][36] . In addition to the DOS, the conductivity has also been studied across this transition 25,30,35,37 and has led to estimates of z and ν for a single Weyl cone 32 , which are consistent with the exponents obtained from DOS calculations. In all of these numerical calculations the CCFS inequality ν ≥ 2/d is well-satisfied. It is important to mention that totally independent from this SM to DM transition, at a much larger disorder strength the Anderson localization transition has been established in some of these models 28,35 . The current work is entirely in the lowdisorder regime (where the SM-DM avoided QCP resides) and has nothing to do with the high-disorder Anderson localization transition from a DM phase to an Anderson insulator phase 28 , which occurs at roughly W l /t = 3.75 for the model under consideration with Gaussian disorder (see the Appendix).
Despite all of this evidence for a stable SM phase and a SM-DM QCP in the presence of disorder, the effects of rare regions call their existence into doubt (and also raise the important and relevant question of why the extensive previous numerical work on the problem always indicates the existence of such a SM-DM QCP). As shown in Ref. 27 through a Lifshitz tail 43 type analysis for the DOS 44 , rare quasi-localized eigenstates ("rare regions") will contribute an exponentially small (in disorder strength) DOS at zero energy, thus making the lowest energy excitations diffusive for arbitrarily weak disorder. Therefore, in the strictest sense, there cannot be a disorder-driven SM-DM QCP in this problem since "both" phases must have nonzero DOS at zero energy although the rare region-induced DOS, being exponentially small, may very well be extremely difficult to discern (or more precisely, the SM-DM transition cannot have the DOS being zero in one phase and nonzero in the other phase as one of its features). These rare eigenstates in the Dirac-Weyl case are distinct from traditional Lifshitz tail states in a band gap (e.g. of a regular semiconductor), as they are only quasi-localized (in contrast to the exponential nature of the disorder-induced Lifshitz band tail states in the semiconductor band gap), with the eigenfunctions falling off at short distances as a power law ∼ 1/r 2 of the distance r from the local extremes of the random potential, and presumably being extended and weakly diffusive at much longer length scales. However, none of the previous numerical studies on the SM-DM QCP has ever observed any signs of these "elusive" rare eigenstates, apart from possibly a large conductance tail in the data of Ref. 32, whose relation to rare events has not been made clear (and which may very well be a finite size effect because finite size systems always have finite conductance). Overall, the numerical data in recent papers seem consistent with the existence of a disorderdriven SM-DM QCP, with a notable agreement between analytical and numerical calculations of the dynamic exponent z(≈ 1.5 within error bars).
In this paper, focusing on a particular lattice model of Dirac (and time-reversal symmetric Weyl) fermions in the presence of short range potential disorder, we address these issues (i.e. both the QCP and rare regions on the same footing) by first finding the rare eigenstates in the SM regime, and then exploring the behavior of the model in the vicinity of the SM-DM avoided QCP. We choose a relatively simple model that has been shown 28,33 to exhibit a sharp SM to DM transition (or crossover), without the additional complications of mass terms. Using (separately) Lanczos 45,46 and the kernel polynomial method (KPM) 47 we provide definitive numerical evidence for the existence of two distinct types of low-|E| eigenstates in the three-dimensional undoped (i.e. Fermi level at the band touching point taken to be the energy zero) Dirac-Weyl systems for weak disorder strengths. Focusing first on the distribution of the first few low-|E| eigenstates, we show for weak disorder that the DOS is well described by "Dirac peaks" (the clean eigenstates that have moved and broadened in energy due to disorder) and an orders of magnitude smaller (i.e. rarer) "background" that fills in between these finite-size Dirac peaks giving a nonzero contribution to ρ(E = 0). We are able to systematically establish that the eigenstates that make up the peaks are perturbatively dressed Dirac eigenstates and the smaller background DOS comes from quasi-localized rare eigenstates. As we show, the peak eigenstates are well described by perturbation theory and are Dirac plane waves weakly distorted by the random disorder potential; they disperse linearly from E = 0. The rare eigenstates are quasi-localized (i.e. the wavefunctions fall off algebraically rather than exponentially) and thus weakly dispersive. Our numerical results indicate that these rare eigenstates arising from disorder (with no clean system analogs) are power-law localized like ∼ 1/r x at short distances r from the site/cluster with the largest disorder strength with the power law x in the range 1.5 − 2.0, in excellent agreement with the analytic prediction (∼ 1/r 2 ). 27 We estimate the zero energy DOS from the background rare region contribution using separately Lanczos and the KPM, finding good agreement between the two methods. Our reason for using two completely independent numerical methods in identifying and quantifying rare region contributions to the DOS is to ensure the accuracy and consistency of our results, given the significance of our findings. Over a range of about four orders of magnitude in the DOS, it well satisfies the rare region form ρ(E = 0) ∼ exp(−a/W 2 ), where W is the amplitude of the random potential. As the disorder strength W increases, eventually there is a crossover to the avoided quantum critical (AQC) regime, where there no longer is a clear separation of the eigenstates between dispersive Dirac states and quasi-localized rare resonances, as the magnitudes of the DOS contributions from the dressed Dirac states and the rare regions start overlapping. In this crossover AQC regime, the DOS far enough away from E = 0 does show a quantum critical form ρ(E) ∼ |E| (d/z)−1 with a z ∼ = 1.5 (Ref. 33), but this scaling behavior is cut off at lower energies. Thus, we conclude that for the model under consideration (and other models with similar symmetry considerations) the SM-DM QCP is converted by nonperturbative effects into an avoided QCP, although the crossover effects of the AQCP manifest themselves at nonzero energies in spite of the QCP itself being suppressed. Our results, taken together with previous work, are consistent with a QCP that is "hidden" by effects that are nonperturbative in the disorder. But a quantum critical regime still exists over a range of nonzero energies, where the rare region correction to quantum critical scaling is small and the nonzero energy behavior of the avoided QCP can therefore manifest itself. The actual size of this crossover region in the energy/disorder space depends crucially on the nonuniversal details of the problem. In other words, the nonzero value of ρ(E = 0) due to the rare eigenstates cuts off the divergence of the correlation length at some length scale ξ RR , thus for length scales ξ < ξ RR , and over the corresponding energy scales, the model looks critical. This is why the previous numerical studies 24,28,32-36 observed an "apparent" SM-DM QCP. Our work thus not only establishes the nonexistence of the disorder-driven SM-DM QCP at the Dirac point due to rare region effects, but also reconciles the large body of existing numerical work finding the existence of such a QCP by showing that the QCP becomes avoided at some large enough length scale (i.e. the correlation length never diverges in the thermodynamic limit) whereas at length scales smaller than this rare region induced cut off length scale, the observed behavior is consistent with a QCP. Fig. 1 gives an overview of the behavior: At weak disorder in the SM regime (e.g. W = 0.4t, where t is the usual nearest-neighbor kinetic hopping amplitude as in Eq. (1)), the DOS is close to the expected ρ(E) ∼ E 2 .
But actually there is a very small nonzero DOS at E = 0 that can not be seen on this linear plot. In the avoided quantum critical regime near W = 0.75t, over a significant range of |E| the DOS is closer to the expected QC behavior of ρ(E) ∼ |E|, although this singularity is always rounded out near E = 0 due to the rare region induced contribution which cuts off the quantum criticality. We are able to quantify how rounded out the singularity is by fitting the low-energy DOS to an analytic form. Then in the DM regime at even larger values of W the nonzero ρ(E = 0) becomes large.
The numerical work presented here establishes three distinct aspects of the disordered Dirac spectra: (1) well below the putative SM-DM transition in the weak disorder regime, there is a nonzero DOS at zero energy; (2) this nonzero DOS arises from the rare regions and is not due to dispersive Dirac quasiparticles, and obeys the expected rare region phenomenology (power-law quasilocalized eigenstates, and exponentially small DOS); (3) this converts the phase transition in to an avoided QCP, which still exhibits a quantum critical regime, but at the lowest energies and the longest length scales becomes an apparently nonsingular crossover between the diffusive metal regime and the regime of a semimetal with these rare quasi-localized eigenstates.
In order to get only one of the two degenerate eigenvalues associated with the conservation of axial charge we construct a two-component model defined as where χ r is a two-component Pauli spinor, the σ µ are the Pauli operators, and there is now only a degeneracy due to time reversal symmetry and the model represents a Weyl Hamiltonian. In the following we will only work with H W and from this point on refer to it as H. We can remove time reversal symmetry by putting twisted boundary conditions on our samples of size L × L × L, so t µ = t exp(iθ µ /L) with −π < θ µ < π. Then there are generally no degeneracies in a finite-size system, and the effect of the random potential at first order in perturbation theory is to simply rigidly move the energies of all these plane-wave eigenstates by the average value of the random potential, which is of order L −3/2 . To remove this leading order finite-size effect, we shift the random potential to always have mean zero: The unshifted random potential V (r) at each site is chosen independently from a Gaussian distribution with zero mean and standard deviation W . We useṼ (r) to denote the shifted random potential with mean zero: We always use energy units with t = 1.
The clean system with W = 0 has a spectrum that consists of discrete levels in any finite-size system. Once we then average over the random potential at nonzero W in the semimetal regime, these discrete Dirac energy levels each give a broadened peak in the DOS. We want to minimize this broadening as much as possible in order to be able to see the rare quasi-localized eigenstates at low energies in between these Dirac peaks. This is the motivation for shifting the random potential to always have zero average, see Fig. 2. This does not change the system at all in the limit of large L, but it changes the finite-size effects on the disorder-averaged DOS, making it easier to clearly see the rare region contributions in spite of their small values. 1t as a function of LE using the shifted random potentialṼ (r) (a) without a twist and (b) with a twist of θx = π/4. The Dirac states that exist in the clean limit are broadened by the disorder, but remain well-separated. Without a twist there is a Dirac peak at E = 0. Applying a twist splits this peak and pushes it away from zero energy, so that there are no states near E = 0. We have checked that in between the peaks there are no states and the flat background seen here is solely an artifact of the KPM. We have also checked (not shown) that this KPM background is independent of NC , provided NC is not too small.
One of the main results of this paper is to directly detect the nonzero DOS at zero energy for weak disorder in the semimetal regime, arising from rare quasi-localized eigenstates. To do this, we use twisted boundary conditions such that in the clean system (W = 0) the DOS of a finite system does indeed strictly vanish at E = 0. Standard periodic boundary conditions for this system unfortunately put Dirac states right at E = 0, thus obscuring this question. Introducing disorder broadens these disorder-averaged Dirac peaks, but for standard periodic boundary conditions the peak remains centered at zero energy, as shown in Fig. 3(a). This connection of a clean Dirac state to its weakly disordered counterpart is made concrete in section III, where we firmly establish that each eigenstate in the peak does represent a (perturbatively dressed) dispersive Dirac state. We can push all of the Dirac states away from zero energy by using twisted boundary conditions (see Fig. 3(b)), which is achieved by t µ = t exp(iθ µ /L) and using periodic boundary conditions. For example, consider a twist of θ = (π/4, 0, 0) as in Fig. 3(b), this pushes the lowest energy eigenstates out to E = ±t| sin(π/4L)| (for W = 0), with no state closer to E = 0.
For the Lanczos calculations that follow we consider odd L and use a twist of θ = (θ x , 0, 0), usually with 0 < θ x < π/2. This is enough to lift the degeneracy of the eigenstates with the four lowest |E|. Focusing on nondegenerate states is preferable for Lanczos, since it has difficulties with degenerate eigenvalues. When we want to push all of the Dirac states as far away from zero energy as possible we use a twist of θ = (π, π, π) with even L, which places the lowest energy eigenstates at ±t √ 3| sin(π/L)|; we find this to be helpful when we use the KPM to estimate the rare eigenstate contribution to the zero-energy DOS. Finally, when we want to estimate ρ(E) while minimizing the finite-size effects at all E, we use KPM and average over all possible twisted boundary conditions, as in Fig. 1.
We study the low energy eigenstates of H using Lanczos and separately the KPM. Lanczos provides accurate estimates of eigenstates and eigenenergies provided the spectrum has no near-degeneracies. Therefore, we focus on the two-component model, usually with odd L and a twist θ x = π/3. We cannot use Lanczos at the largest disorder strengths, due to the spectrum becoming too dense near zero energy so Lanczos misses states due to the near-degeneracies. In the semimetal regime, the four lowest-energy Dirac states for this twist and odd L are at energies near ±E 0 and ±2E 0 , where E 0 ∼ 1/L depends on both L and W . Even though disorder breaks the particle-hole symmetry (E → −E), for weak disorder it is only weakly broken and Lanczos on H 2 combines states at the (approximately) same |E| in to the same peak in ρ(|E|), as shown in Fig. 4. In order to space out these eigenvalues we also do Lanczos on (H − (E 0 (W, L)/4)) 2 . This puts the first four states near effectively separating the peaks for each Dirac state, provided the width of each peak is not too broad.
Focusing on N u eigenstates from Lanczos the low energy average DOS can be computed from their distribution, for N R disorder realizations, where E i (r) is the ith eigenvalue of the rth disorder realization. We often get the  . The two-component nature of the DOS, consisting of Dirac peaks and a smooth "background" is clear over an intermediate range of disorder W in the semimetal regime. The background DOS extending down to zero energy is detected for W ≥ 0.5t. As W increases, the Dirac peaks broaden and eventually the clear distinction between peaks and background is lost at larger W in the avoided quantum critical regime.
N u = 4 lowest states. A few comments about the definition of ρ low (|E|) are in order: First, this definition implies that the DOS is normalized as the number of states per volume per dE, and therefore will have the same meaning as the full DOS for a particular energy E. The (full) average DOS is defined as where D = 2L 3 is total number of states for the twocomponent model. Second, the low energy DOS in Eq.
(3) is only an accurate estimate of the full DOS for energies |E| ≤ E * Nu where E * Nu = min {r} |E Nu (r)|, i.e. the minimum value of the largest (N u ) eigenvalue Lanczos has computed. As a result, whenever we show ρ low we will also plot vertical dashed lines to mark E * Nu (W, L) (apart from Fig. 2 where this is not an issue due to the very weak disorder). For |E| > E * Nu Lanczos begins to miss some energy eigenvalues and this low energy estimate of the DOS will be depleted relative to the full ρ(E). Lastly, since Lanczos has individual eigenvalue resolution (as opposed to the KPM), we bin the results to generate a smoother estimate of the DOS.
For the KPM calculations presented in section IV we consider twisted boundary conditions with θ = (π, π, π). The technical details of the KPM can be found in Ref. 47. KPM essentially trades off computing eigenvalues of H for directly computing ρ(E) via an expansion in terms of Chebyshev polynomials to an order N C , and uses a kernel to filter out Gibbs' oscillations due to truncating the expansion. (Avoiding direct diagonalization to calculate the eigenenergies allows the KPM to go to very large system sizes to directly compute the DOS, which would be inconceivable within an exact diagonalization.) Here we use the Jackson kernel 47 , which amounts to replacing the delta function in Eq. (5) with a normalized Gaussian with a standard deviation σ = π/N C , (in units of the bandwidth ∼ 2( √ 3t + W )) so that the DOS over the full bandwidth remains normalized to ρ(E)dE = (total number of states per volume). In the calculations that follow we are using N C = 1024 unless otherwise stated. When we average over samples, we use the same energy grid for all samples. In order to effectively use the KPM to study the rare region contribution to the DOS, we find it essential to consider two issues: First, due to this artificial broadening, using the twisted boundary conditions to push the Dirac peaks as far away from zero energy as possible is helpful, so that they do not contaminate the estimate of ρ(E = 0). Second, even if there is a strictly zero DOS at a particular energy E, the KPM will give a non-zero number for ρ(E) and will thus give an artificial "background" to the KPM DOS, for example see Fig. 3. As a result of this artificial background, even in the presence of the twist we find that the KPM cannot accurately determine ρ(0) for the smallest W of interest. Therefore, we can use the Lanczos to obtain ρ(0) for weak disorder and the KPM for large disorder, whereas for intermediate disorder strengths we find that the estimates from the two methods do match consistently, which is an important numerical check for our results.

III. EIGENSTATES OF H
In this section we study the nature of the eigenstates of H. As we discuss below, in the semimetal regime at weak disorder we find two qualitatively distinct types of eigenstates that give separate contributions to the DOS for finite samples. In this regime, we find that the DOS can be separated into "peaks" and a "background" that lies in between the peaks. This separation is useful as it will allow us to study the eigenstates that make up each contribution separately. This is shown clearly in Fig. 4 for L = 25 and 10,000 disorder realizations using Lanczos on H 2 for the first four lowest energy eigenstates. Note that H 2 has put the Dirac peaks that are at the (approximately) same value of |E| on top of each other, so in this figure the four states at E ∼ = ±E 0 , ±2E 0 produce two peaks.
For very weak disorder we see the two expected Dirac peaks with a "background" DOS developing between the peaks. For larger disorder strengths in the semimetal regime the Dirac peaks remain and in addition we find the background is detected all the way down to zero energy. We expect that the low energy tail is orders of magnitude too small for W/t = 0.30 and 0.40 to be observed for these system sizes and this number of disorder realizations, but is still actually present at any nonzero W . We find this background DOS is an increasing function of |E|. For still larger disorder in the AQC regime the distinction between "peaks" and "background" is eventually lost. For these larger disorder strengths we expect that the excitations are all diffusive and there are no longer well defined dispersive Dirac excitations; this occurs where the Dirac peaks are no longer visible (W ∼ = 0.7t). For lower disorder strengths, the low energy background DOS represents quite rare eigenstates, e.g. for W = 0.5t the magnitude of the DOS at the peak versus the low energy (rare region) background is separated by almost four orders of magnitude. This explains why earlier numerical work invariably failed to find any rare region induced background DOS, thus concluding (erroneously) that the system remains a semimetal with zero DOS at E = 0 up to the critical disorder strength.
Previous numerical studies of the SM-DM transition 24,28,33-36 used periodic boundary conditions and even L, which produces a very strong finite-size effect on the zero-energy DOS in the semimetal. In the current model with Dirac points occurring at momenta commensurate with the lattice, even L and periodic boundary conditions place the lowest energy eigenstates into the Dirac-peak centered around E = 0. In models where the Dirac/Weyl cone is located at momentum incommensurate with the lattice, this is not the case. However, without carefully choosing L or using a twist, Dirac/Weyl states will inevitably come sufficiently close to zero energy producing a strong finite size effect at E = 0 (this is straightforward to see from our results since it is essentially equivalent to considering a small twist in our current model). By pushing the Dirac states away from zero energy we have now allowed the background DOS at zero energy to be visible in the semimetal regime. In the remainder of this section we will now study the two different types of eigenstates separately, i.e. the rare eigenstates that contribute to the low energy background DOS and the typical Dirac states that make up the peaks in the DOS. In the next section we will estimate the background contribution to the zero energy DOS.

A. Rare eigenstates
In Ref. 27 the theory of rare quasi-localized eigenstates in three-dimensional Dirac systems has been derived from a Lifshitz tail type formalism with the essential idea that the rare region effects in the Dirac-Weyl gapless systems are basically like the resonant versions of Lifshitz tail states which reside in energy band gaps. It was found that the eigenstate that corresponds to a rare event, i.e. a disorder configuration that has either a site or a small group of sites that has a very large disorder strength, is a quasi-localized resonance that decays from the site r max with maximal disorder in a power law fashion like ψ(r) ∼ 1/r 2 for r/L ≪ 1 and r ≡ |r−r max | more than one lattice spacing. It is important to stress that the existence of these rare eigenstates of H is non-perturbative in the disorder strength, and hence is outside the scope of the perturbative field theoretic analysis of the SM-DM QCP discussed in the literature.
We will now study the properties of eigenstates in a particular rare disorder realization that gives rise to states in the background DOS. In this subsection, we focus on a sample of size L = 25 and a disorder strength W = 0.5t. By varying the twist θ x in the x-direction we can determine the dispersion of an eigenstate in the "mini" Brillioun zone for momentum −π/2 < θ x < 3π/2. Focusing on the lowest four eigenstates of H 2 , we determine the sign of each eigenvalue of H from its corresponding eigenvector, and construct the dispersion of these four eigenstates in both positive and negative energies as shown in Fig. 5. There are two weakly dispersive and thus quasi-localized states and two dispersive Dirac states, and these states hybridize near avoided level crossings. The states come in pairs with opposite spin. For twist θ x = 0 and π the system has time reversal symmetry and thus degenerate Kramers doublets. In this sample, the rare states have a small negative energy, but (color online) Projected probability density z |ψ(x, y, z)| 2 versus x and y for a weakly-dispersing rare state with L = 25, W = 0.5t, and a twist θx = π/2. Note the system has periodic boundary conditions (with tx = t exp(iθx/L)), and we have set the lattice spacing to unity.
among samples with such states, the energy is smoothly distributed through zero energy, resulting in a nonzero contribution to the zero-energy density of states. Now that we have determined how this rare state disperses, we turn to the magnitude of the wavefunction ψ(r) ≡ |ψ a (r)| 2 + |ψ b (r)| 2 where a and b label the two spinor components for the lowest-|E| eigenstate with a twist that makes the state the most weakly dispersing, i.e. at θ x = π/2. For plotting purposes we show z |ψ(x, y, z)| 2 , which clearly shows a (quasi) localized wave-function, see Fig. 6. We find the location in real space where two neighboring sites have a large disorder strength V i ∼ 3W at the same location r max where the wave-function's magnitude is maximal. The probability of this disorder configuration is quite rare, relative to the probability of a typical configuration (V i ∼ W ) it occurs with a probability ∼ exp(−9), and therefore this is indeed a rare eigenstate.
We define the decay of the wave-function from its maximal value by computing ψ(r) ≡ ψ(|r − r max |) (for |r µ − r µ max | < L/2 respecting the periodic boundary conditions). In Fig. 7(a) we show the scatter plot of the decay of the wavefunction from its maximal value, which indicates a power law trend in the data. We then discretize the r axis into bins and then average the value of ψ(r)/ψ(r max ) in each bin, which yields ψ bin (r) as shown in Fig. 7(b) for W = 0.5t. To demonstrate the sample to sample variations we also show ψ bin (r) for W = 0.6t in Fig. 8 for two different disorder realizations that give rise to distinct quasi localized eigenstates. We now reach one of our main results, where over a range of r we find the (binned) rare wave-function decays like with a power law exponent x that varies realization to realization of disorder in the range 1.5 − 2.0, which is in excellent agreement (within numerical accuracy) with the analytic prediction of ψ(r) ∼ 1/r 2 in Ref. 27. It is interesting to contrast these eigenstates with states in the Lifshitz tail of the DOS in the presence of a band gap. The latter are exponentially localized 48 at a site with a very large disorder and contribute an exponential tail to the DOS near the band edge 43 . Here, there is no band gap, and as a consequence these rare eigenstates are only power-law bound on these short length scales, and it is in this sense that they are only quasi-localized. Intuitively, these rare states "pull" some spectral weight out of the Dirac bands and can place that weight at arbitrarily low energy, thus contributing to a background DOS that remains nonzero through E = 0. Thus, the rare regions destroy the simple distinction between the SM phase (with zero DOS at zero energy) and the DM phase (with nonzero DOS at zero energy) as the DOS is always nonzero albeit very small for weak disorder. The system is unstable to any disorder which immediately produces a nonzero DOS at zero energy. (color online) Projected probability density z |ψ(x, y, z)| 2 versus x and y for a bi-quasi-localized rare state in the low energy tail of the DOS with L = 25, W = 0.66t, and a twist θx = 0.325π. Note the system has periodic boundary conditions (with tx = t exp(iθx/L)), and we have set the lattice spacing to unity.
For increasing disorder strengths that remain in the SM regime, the probability to generate more then one rare region increases, which makes it increasingly likely to find multiple quasi localized power law states in a single wavefunction per sample. Again focusing on an eigenstate that contributes to the low energy background DOS for W = 0.66t, in Fig. 9 we show z |ψ(x, y, z)| 2 which clearly reveals the existence of a bi-quasi-localized wave function. Since each wavefunciton falls off (roughly) as 1/r 2 from the sites (r 1 and r 2 ), the overlap of these two quasi-localized peaks produces a non-zero tunneling matrix element that goes as t RR (r 1 − r 2 ) ∼ 1/|r 1 − r 2 | 2 . Therefore it is natural to expect this tunneling will produce a diffusive metal where the conductance is mediated by hopping between these rare regions of large probability amplitude 27 . We do expect that such wave functions are also generated at much weaker disorder, however their probability is so small that they are essentially never found in these size samples.

B. Perturbatively dressed Dirac eigenstates
We now consider the eigenstates that make up the low energy Dirac peaks in the DOS. For W = 0.3t and L = 25, by varying the twist we again determine how the four lowest-|E| states disperse in the mini-zone for a typical sample. As shown in Fig. 10, we find that the states disperse linearly, just as they do in the absence of   Fig. 4 into two) as a function of the disorder strength W , the dashed line is a fit to the perturbative form bW 2 with fit parameter b. We find the fit works reasonably well for W/t ≤ 0.1 and this FWHM increases more strongly than quadratically in W for larger disorder strengths.
the random potential. The only visible effects of the random potential are a small renormalization of the Fermi velocity that we discuss below, and weakly-avoided level crossings at θ x = ±π/2, where the resulting eigenstates are standing waves. The probability density in one of these standing wave eigenstates is shown in Fig. 11. Here we can see that the eigenstate is very regular, with only weak "noise", due to the perturbatively irrelevant disorder. We now focus on the quantitative properties of the lowest-|E| Dirac peak at twist θ x = π/3, and their dependence on W and L (see Appendix A for the perturbative analysis which is consistent with our numerical results for states in the Dirac peaks). By restricting the random potential to have sum zero, we have eliminated the first-order-in-W perturbative effect. At order W 2 , there is level repulsion from all other momenta, which is dominated by the many states that are far away in energy, since the DOS is so small at low |E|. The net effect of all this level repulsion is to reduce the Fermi velocity at order W 2 , because positive (negative) energy states have stronger repulsion from the other positive (negative) energy states, since they are closer in energy, and thus the mean energy is pushed down (up) by the level repulsion. This is illustrated in Fig. 12(a), where we see a ∼ W 2 suppression of the energy fits well over most of the semimetal regime. The random component of the level repulsion gives the sample-averaged Dirac peaks a linewidth ∼ W 2 /L 2 for small W (see Fig. 12(b) for the full-width-at-half-maximum (FWHM) Γ versus W ).
Here we find the dependence on W fits this quadratic behavior only for quite small W , and the linewidth increases faster than ∼ W 2 throughout most of the SM regime.
We now turn to the system size dependence of the |E|peaks. For the linearly dispersing Dirac excitations we find that each peak's energy follows the leading 1/L dependence inherited from its clean limit behavior as shown in Fig. 13(a), with an exponent x from E ∼ 1/L x that varies from 0.995 to 1.010 (see the inset of Fig. 13(a)). We find this behavior clearly up to a disorder strength of about W ∼ 0.62t (not shown), beyond which cleanly identifying the energy of the peak above the (relatively) large background is no longer possible. We have also checked that this holds for the second peak, as well, up to W ∼ 0.62t (not shown). The FWHM Γ of the first peak is shown in Fig. 13(b), which is well described by the perturbative result Γ ∼ 1/L 2 only for very weak disorder strengths W ≤ 0.1t, which is consistent with where Γ deviates from the perturbative expectation (∼ W 2 ) as in Fig. 12(b). For W > 0.1t we find a systematic decrease of the exponent x in Γ ∼ 1/L x .
Thus we have shown that in the semimetallic regime, most properties of the eigenstates that make up the lowenergy Dirac peaks in the DOS are well described by treating the Dirac eigenstates perturbatively in the disorder strength. These Dirac peaks in the spectrum survive up to a disorder strength W ≈ 0.6 t. The one property that is not well captured perturbatively is the dependence of the width of the Dirac peaks on W , which grows faster than ∼ W 2 in most of the SM regime.
For larger disorder strengths the model enters the avoided quantum critical regime, with a substantially nonzero ρ(0), and there is no longer a clear distinction between the peak and the background contributions to the DOS. In this regime the rare eigenstates (that make up the background) are no longer rare at all and become typical eigenstates, as the background fills in between the peaks. It is both compelling and consistent that in our microscopic study of these Dirac peaks we find that they are no longer clearly part of the excitation spec- trum for disorder strengths W > 0.62t, where the earlier KPM study using periodic boundary conditions on this model found the diffusive metal regime. Thus, the SM-DM crossover behavior survives the existence of the rare regions although the SM-DM QCP itself is destroyed by the rare regions!

IV. DENSITY OF STATES
We are now in a position to estimate the rare eigenstate contribution to the zero-energy background DOS. It is important to emphasize that this is in contrast to the estimate of the Dirac peaks' contribution to ρ(0) in Refs. 28  , and W = 0.58t (c). The vertical dashed lines mark E * Nu where for larger energies the DOS may be underestimated due to this calculation not getting the fifth and higher states. The label for each L is the same for each plot. We find a clear L-independent low energy tail, which we take as an estimate of ρ(0). We also find that the Dirac peaks remain and continue to sharpen up for increasing L in this range of L and W , as detailed in Figs. 12 and 13. (H − E 0 (W, L)/4) 2 with a twist θ x = π/3 and 10,000 disorder realizations for each value of W and L. We then move onto the results using the KPM on H with a twist of θ = (π, π, π) and 1,000 disorder realizations. The estimates of the low energy density of states from Lanczos on (H − E 0 /4) 2 are shown in Fig. 14 for disorder strengths W/t = 0.50, 0.54, 0.58 which are all in the semi-metal regime. We find that the background contribution to the DOS develops an L-independent low-energy tail, which is one of our main results. Therefore this estimate of ρ(0) is nonzero in the thermodynamic limit, albeit quite small for a weak disorder strength. We also find that in this regime the Dirac peaks remain and continue to sharpen for increasing L; the energies and widths of them are shown in Figs. 12 and 13 and discussed in section III. The orders of magnitude difference between the background and the Dirac-peak contributions to the DOS indicates the difficulty in trying to observe the rare region contribution to the DOS at zero energy without appropriately choosing the boundary conditions to explicitly separate these distinct contributions to the DOS. We should note that these Lanczos estimates of the low energy DOS are actually of ρ(E 0 (L)/4), thus not at strictly zero en-ergy. But we find that in this very low energy range the DOS has only a weak energy dependence, so this is not a significant difference, especially compared to the roughly four order of magnitude range of the DOS as we vary W .
When estimating ρ(0) using KPM, we would like to push the Dirac peaks as far away from zero energy as possible. The Dirac peaks are broadened both due to disorder and due to using finite N C in the KPM, and we want to minimize any contribution from this broadening to the estimated ρ(0). This is achieved by using a twist of θ = (π, π, π) in each direction and even L. It is important to also remember that the KPM introduces an artificial KPM background (as shown in Fig. 3) which at small W contaminates our estimate of the true background DOS. Therefore we cannot extend the KPM estimates of ρ(0) to as small W as we have for Lanczos. As shown in Fig. 15, using the KPM for disorder strengths W/t = 0.600, 0.625, 0.650, we find a flat low energy background contribution to the DOS that is L-independent and extends all the way to E = 0. Similar to the Lanczos data we still observe the Dirac peaks separating the smooth DOS at higher energy from the flat background, although for W > 0.6t these peaks are being rounded out in to "shoulders". For energies above these Dirac peak/shoulders we find that in this AQC regime the DOS is close to the quantum critical (QC) form ρ(E) ∼ |E|, which is in agreement with the data in the absence of a twist 33 . Thus, the crossover effects of the QCP survive, but the QCP itself does not.
The Dirac peaks in the SM regime are separated by 1/L (see inset of Fig. 15(b)) and are a finite size effect, whereas the DOS is converged in L both at high energy and near zero energy. We can remove this finite size effect by averaging over the twist, which we do by generating a random twist vector θ = (θ x , θ y , θ z ) where each θ i is a random twist for each disorder realization, uniformly distributed between 0 and π; thus we average over the twist and disorder. This is displayed in Fig. 16 (a) for 1, 000 disorder/twist realizations and a disorder strength W/t = 0.6. For system sizes L ≥ 37 we find the data is well converged in L for all energies, and as a result the data clearly displays three regimes in energy: the DM regime at the lowest energies, where the DOS does not depend on E; the SM regime at intermediate energies, where ρ(E) ∼ E 2 ; and the QC regime at higher energy, where roughly ρ(E) ∼ |E|. As shown in Fig. 16 (b): At smaller W , the QC regime disappears by moving up to near the cutoff. At larger W , the size (in energy) of the DM regime grows and the SM regime disappears, leaving a direct DM to QC crossover, and at much higher W the QC regime again disappears by moving up to near the cutoff. It is important to emphasize that the DOS is always smooth through E = 0 and L-independent for our largest values of L. We can directly characterize this by expanding the DOS as ρ(E, W ) ≈ ρ(0, W ) + a(W )|E| 2 + . . . . As shown in Fig. 16(c), a(W ) rises smoothly developing a finite peak, with no divergence and no sign of any QC singularity at E = 0. The location of the peak value of a(W ) provides an estimate of the avoided QCP: we find W c ≈ 0.75t, which is consistent with the cross over regime in energy in Fig. 16 (b). The implications of this are twofold: First, estimating the location of the QCP based on an apparent vanishing of ρ(0) actually underestimates W c because (as we have shown) ρ(0) is always non-zero. Second, there is no indication at all of any singularity in the DOS in this system at any critical value of W and the QCP is clearly rounded out ("avoided"). This establishes that our numerical results are inconsistent with the DOS being expressed as a sum of two independent terms, a singular one arising from the quantum criticality and a smooth background contribution from rare regions. This overall behavior is illustrated in the schematic disorder versus energy crossover diagram in Fig. 16(d).
Using the two estimates we now have from the Lindependent part of the background DOS we fit our data to the rare-region form which fits remarkably well over four orders of magnitude of ρ(0) going down to W = 0.48t from Lanczos and up to W = 0.8t for the KPM, with again no sign of any singular behavior at a critical value of W , as shown in Fig. 17(a). We find good agreement between Lanczos and KPM estimates of ρ(0) which implies that the background DOS is twist independent. We find that the slopes of the fitted lines in Fig. 17 (a) are in good agreement and match to within ∼ 1% with a small offset between the two (about a 10% difference). We attribute the systematic underestimate of ρ(0) from Lanczos due to this method missing nearly degenerate eigenvalues that can arise from squaring the Hamiltonian. In addition, for small disorder strengths we see the KPM estimate starts to "peel off" systematically from the fitted line, we attribute this to the artificial KPM background in the DOS setting a lower bound to the value of ρ(0) that we can accurately estimate with the KPM.
To conclude this section, we now discuss the various estimates of ρ(0) using Lanczos, the KPM with no twist, and the KPM with a twist of π in all directions, as shown in Fig. 17 (b). Deep in the DM regime at high W , we expect that ρ(0) should be completely independent of the twist and L (note this region cannot be reached with Lanczos as the spectrum is too dense), which is in good agreement with the numerics shown in Fig. 17 (b). For weaker disorder strengths the DOS from the KPM with a twist and the Lanczos estimates exponentially decrease and remain L-independent, becoming orders of magnitude smaller than the L-dependent value of ρ(0) without a twist. As we illustrated in Fig. 3(a), this is because a Dirac peak sits at E = 0 if there is no twist, making this a very strong finite-size effect for ρ(0) in the semimetal regime. The onset of this strong finite-size effect occurs near the avoided quantum critical point. The avoided quantum critical point also affects the DOS away from zero energy, giving it a quantum critical scaling of roughly ρ(E) ∼ |E| over a range of energy. A crossover diagram summarizing this is shown in Fig. 16 (d). We emphasize that although this crossover diagram is schematic, all aspects of it are obtained from our numerical results presented in this work.

V. DISCUSSION AND CONCLUSION
Our results have demonstrated that in the semimetal regime the quasi-localized rare eigenstates live on a continuum of low energies that contribute a nonzero L independent low energy DOS. These rare eigenstates do not live in isolation and in principle there can actually be several per sample (here we have demonstrated a pair of these resulting in a bi-quasi-localized wavefunction as shown in Fig. 9) with a non-zero tunneling matrix element that falls off with separation between the peaks (r) like ∼ 1/r 2 . All of these results are suggestive that these rare states are not fully localized. Therefore we now turn to the low energy level statistics to directly address this Avoided QCP Note that the value of ρ(0) here is larger than for a twist of θ = (π, π, π), as averaging over the twist does mix together the Dirac and rare states contributions to the zero energy DOS. (c) Fit parameter a(W ) as a function of W extracted from fitting ρ(E) − ρ(0) to a(W )|E| 2 in the low energy limit from the KPM averaged over the twist and disorder for 1, 000 realizations with a linear system size L = 71. We find a(W ) is a smooth function of W and provides an accurate estimate of the avoided QCP Wc = 0.75t, thus we find the DOS is not described by the combination of a smooth background and a critical part the critical point is sufficiently rounded out. (d) Schematic crossover diagram as a function of energy and disorder strength. Despite the existence of rare region effects we still find semimetal and quantum critical regimes exist, albeit at nonzero energies. The quantum critical scaling regime is "anchored" by the avoided QCP and consistent with the perturbative RG analysis.
question. We compute the adjacent gap ratio defined as where δ i = E i+1 − E i is the level spacing between neighboring energy eigenvalues. Here we focus on the center of the band and only compute r i for the lowest four |E| eigenvalues through exact diagonalization on even L samples with periodic boundary conditions, so that the level spacing will capture a mixture of both Dirac and rare eigenstates. Focusing on weak disorder where we can find the background DOS, as shown in Fig. 18 we find that the disordered average r ≈ 0.53 (to within numerical accuracy) for the biggest L we have considered and therefore the level statistics satisfy the Gaussian orthogonal ensemble (GOE) 49 . We also find that r is unaffected by crossing the avoided QCP at W ≈ 0.75t. This establishes that the low energy eigenstates have a non- zero level repulsion and thus are not localized eigenstates. Therefore we can safely conclude that the quasi-localized eigenstates that fall off in a power law fashion are not localized, in strong contrast to the exponentially localized Lifshitz states (that live in a band gap or band edge).
In this work we have considered the effects of potential disorder on a three-dimensional model that possesses an axial symmetry, which is pertinent to describing the bulk physics in various Dirac semimetals as well as timereversal-invariant Weyl systems. We expect our results to be broadly applicable to models with the same symmetries as H D or H W . In this regard, Ref. 33 found that ρ(0) was numerically independent of tuning either potential, axial, or mass disorder for three-dimensional Dirac fermions, however now our work establishes that this is true with regards to the Dirac eigenstates only (as this was what was being computed in Ref. 33). Nonetheless, since the model with potential or axial disorder has a con- tinuous axial symmetry, they can both be written in the form of H W while for mass disorder they cannot 28,33 and thus our results also apply to axial disorder. It will be interesting to see if this observation remains true with regards to the background (rare eigenstate contribution) to the DOS. It will be exciting to explore disordered models with other symmetries, such as, for example, cases with disorder that preserves particle-hole symmetry. Perhaps there are other models where the semimetal is stable to disorder and a true phase transition out of the semimetal occurs at some nonzero disorder strength. This is, however, well beyond the scope of the current work where we consider the canonical (and by far the most studied) model of short-range potential disorder in the context of disorder driven SM-DM QCP, finding that the QCP does not exist (or becomes avoided) due to rare region effects, but the crossover effects of the QCP may exist in the quantum critical fan region.
Our results raise the question: Is there a field-theoretic description of these rare eigenstates and this avoided critical point? It is now clear that the self-consistent Born approximation and perturbative RG are not capturing the crucial but nonperturbative effects of the rare disorder configurations. It is important to mention that the analytic theory of these rare eigenstates has been obtained following a "Lifshitz tail" analysis and Lifshitz tail eigenstates can also be taken into account within a field theory context and appear as an instanton config-uration 50 . It is therefore suggestive that the rare eigenstates that we have studied here will contribute some sort of non-pertubative instanton configuration that fundamentally changes the perturbative result (i.e. the self-consistent Born approximation) and the associated renormalization group analysis based on loop expansions. Constructing this effective action incorporating the existing field together for both Dirac 33 and Weyl 29,38 fermions should provide an effective theory for the avoided QCP and remains an important open question for the future.
The direct consequences of our findings on the nontrivial topological properties of clean Dirac and Weyl semi-metals such as the surface Fermi arc states 16 and the axial anomaly 51-53 is an interesting and open question. The absence of a bulk gap and the presence of these rare quasi-localized eigenstates may provide a scattering channel from surface arc states into the bulk endowing them with a non-zero quasiparticle lifetime and dissipative transport properties 54 at sufficiently low energies. With regards to the axial anomaly that has been indirectly observed through a measurement of the longitudinal magneto resistance 55-57 , we expect this does survive for reasonable transport time scales 58 as the the quasiparticle lifetime at low energies and weak disorder goes as τ ∼ 1/ρ(0) and will therefore be exponentially large in the strength of disorder. However, our results do establish that any perturbative treatment of the problem that is carried out at nonzero Fermi energy 59 cannot be extended down to E = 0 (i.e. to the Dirac or Weyl cone) because of the existence of non-perturbative rare states. Lastly, the physics of the axial anomaly at larger fields in the quantum limit is unchanged by our findings because this is well described by quasi one dimensional dispersive states that host their own chiral anomaly in one dimension 60 . Both of these questions are sufficiently interesting and the effects of non-perturbative states upon them are sufficiently nuanced that they warrant their own separate study well beyond the scope of the present work.
It is interesting to compare our results for the avoided QCP with that of various strongly-correlated systems (such as heavy fermion metals 61,62 or cuprate superconductors 63 , where the evidence of a QCP is quite striking) where broken symmetry phases set in (such as superconductivity) and mask the zero temperature transition. In these systems the quantum critical features are observed within the quantum critical fan and have a strong effect on finite temperature thermodynamic and transport properties. With this in mind, and the schematic crossover diagram in Fig. 16 (b) that contains a quantum critical fan that is anchored by the avoided QCP, we still expect that if experiments on Dirac and Weyl semimetals can be tuned to the (zero energy) Dirac point thermodynamic signatures of the avoided QCP should show up in the crossover regime, e.g., a specific heat varying like ∼ T 2 . Upon lowering the temperature this power law will eventually cross over to ∼ T due to the rare regions masking the QCP. Thus, in the current problem, the QCP is truly avoided (rather than 'hidden') since there is no way, even as a matter of principle, to think of a situation to restore the QCP since disorder is the tuning parameter both for creating the QCP and for producing the rare regions destroying the QCP. This is conceptually somewhat different from the situation with heavy fermions or cuprates where the origin of the superconductivity might be distinct from the origin of the QCP, at least as a matter of principle, so one can imagine suppressing the superconductivity (e.g. by applying a strong magnetic field) to restore the QCP. In the Dirac-Weyl system, our current work definitively establishes that the disorder-driven SM-DM QCP does not exist as it appears to have been suppressed by the finite density of states contributed by the rare regions. What does survive, however, is the crossover effect of the QCP which should produce effective scaling behavior provided one is at reasonably high energy (i.e. high temperature and/or high frequency). Such an apparent 'effective scaling' behavior, numerically observed in many earlier theoretical studies, has led to the erroneous conclusion on the existence of a disorder-driven QCP in Dirac-Weyl systems, which our current work establishes as being nonexistent since it is avoided at the lowest energy (or the largest length) scale.
To conclude, we have studied the effects of rare eigenstates for Dirac and Weyl semimetals in the presence of potential disorder. Using Lanczos (on H 2 ) and KPM (on H) with twisted boundary conditions we have established a systematic method to isolate and study the effects of rare regions in detail. We have shown that for weak disorder the model under consideration possesses two classes of eigenstates. Consistent with the perturbative irrelevance of the disorder, the first are Dirac eigenstates that are well described by perturbation theory in the random potential. These eigenstates disperse linearly with a wavefunction that is qualitatively consistent with a Dirac plane wave state weakly perturbed by the random potential. The second class of eigenstates are the rare eigenstates that are very weakly dispersive and whose wave function is power-law quasi-localized near a site with strong disorder strength. These eigenstates contribute a background DOS that extends all the way to zero energy and is exponentially small in the disorder strength W . As a result of this non-zero DOS at zero energy, the expected semimetallic regime with DOS ρ(E) ∼ E 2 only exists at energy scales above that set by the rare regions contribution to the DOS, and the apparent semimetal to diffusive metal QCP is avoided, pushing the quantum critical regime to nonzero energy.

VI. ACKNOWLEDGMENTS
We thank Pallab Goswami and Rahul Nandkishore for various discussions and collaborations on related work. We also thank Victor Gurarie, Leo Radzihovsky, and Sergey Syzranov for discussions and for comments on a draft manuscript. This work is partially supported by JQI-NSF-PFC, LPS-MPO-CMTC, and Microsoft Q (JHP and SDS). DAH is the Addie and Harold Broitman Member at I.A.S. JHP acknowledges the hospitality of the Aspen Center for Physics (NSF Grant no. PHY-1066293) where some of this work was completed. We acknowledge the University of Maryland supercomputing resources (http://www.it.umd.edu/hpcc) made available in conducting the research reported in this paper. JHP would like to thank Zhentao Wang for help with the Lanczos package ARPACK.

Appendix A: Perturbation Theory
Here we determine the leading finite size and disorder effects from a perturbative analysis in the strength W of the random potential. For simplicity, we focus on the two-component model (3). We use odd L and include twisted boundary conditions such that the system has no degeneracies at zero disorder, where the eigenenergies and eigenfunctions are E (0) k,± = ±t sin 2 k x + sin 2 k y + sin 2 k z (A1) where the φ ± k are normalized two-component spinors, k 0 is the wavenumber allowed by the twisted boundary conditions that is closest to zero, and l, m, n are integers.

First order correction to the eigenfunctions
We have chosen the random potentialṼ (r) to always sum to zero over all sites, so since there are no degeneracies at zero disorder the first order corrections to the eigenenergies all vanish. The first order contribution to the eigenfunction ψ k,± (r) = ψ (0) k,± (r) + ψ (1) k,± (r) + ... is given by ψ (1) k,± (r) = q =k,s ψ (0) q,s (r) where s is summed over + and −. For nonzero W , this sum is infrared divergent in the limit of large L at energies away from the Dirac point, due to the small energy denominators. This reflects the fact that the mean free path is finite away from the Dirac energy, and those eigenfunctions are highly changed by the scattering when L exceeds the mean free path. However, if we look at the eigenstates with energies closest to the Dirac point, as we do in this paper, then the vanishing of the DOS as ∼ E 2 suppresses this infrared divergence and the sum is instead dominated by typical other states with energies far from E (0) k,± . As a result, the first-order correction to these eigenfunctions is random with only short-range correlations and a relative magnitude that is ∼ W and independent of L at large L.

Second order correction to the eigenenergies
The second order contribution to the eigenenergy is given by We are interested in states near the Dirac point, so let's look in particular at E (2) k0,+ . The contributions from typical momenta q that are far from the Dirac points cancel to leading order: Let's look at the contributions from the 4 other states at momenta q = k 0 ± Q. These two momenta are very close to equal and opposite, and as a consequence the corresponding eigenspinors φ + k0±Q are nearly identical to φ − k0∓Q , the corresponding energy denominators are of opposite sign and of nearly equal magnitude and magnitudes of the matrix elements are nearly identical. Thus the level repulsion from the higher energy states almost exactly cancels that from the lower energy states. What remains from these 4 other states is a random energy shift of order W 2 /L 4 that on average lowers this positive eigenenergy E k0,+ by ∼ W 2 /L 4 . When these are summed over all L 3 other momenta, this gives an average decrease in this eigenenergy by ∼ W 2 /L, which gives a decrease of the Fermi velocity by ∼ W 2 . The random contribution summed over all these typical other states to the eigenenergy is smaller by a factor of L 3/2 , so is ∼ W 2 /L 5/2 . The contribution to E (2) k0,± from other states that are nearby in energy is random and ∼ W 2 /L 2 , so does not contribute to the shift of the Fermi velocity in the limit of large L, but does dominate the random energy shift from typical other states. Thus, when averaged over samples, this "Dirac peak" has linewidth ∼ W 2 /L 2 at second order in W . We only consider a few sites N s ≪ V to improve the statistics and . . . denotes a disorder average. Due to the low DOS in the SM regime the typical DOS is not well suited to study the lack localization of the quasi localized rare states and therefore for this purpose we have used level statistics as shown in the main text in Fig. 18. For large disorder the average DOS is sufficiently large and therefore the typical DOS is well behaved (see Fig. 19 (a)). We use periodic boundary conditions and even L as for these large disorder strengths twisted boundary conditions have a negligible effect. We focus on the localization transition at E = 0, as the mobility edge has been shown to be relatively standard 28 , i.e. it starts at the band edge and decreases in |E| for increasing W .
The results for large disorder and various systems sizes with a KPM expansion order of N C = 4096 are shown in Fig. 19(a), we find ρ t (0) is well converged for L ≥ 30. To study the localization transition we fix the linear system size to L = 40 and vary the KPM expansion order. By extrapolating ρ t (E = 0) to zero in Fig. 19(b) we find an estimate of the localization transition W l as a function of N C , extrapolating to N C → ∞ yields an estimate of the localization transition at W l /t = 3.75 ± 0.25. This places the standard Anderson localization transition for E = 0 at a much larger disorder strength then the avoided QCP.