{'fraction_non_alphanumeric': 0.03409374340476967, 'fraction_numerical': 0.010245390533566125, 'mean_word_length': 4.3018004272200185, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Rare region effects dominate weakly disordered 3D Dirac points

We study three-dimensional Dirac fermions with weak finite-range scalar potential disorder. In the clean system, the density of states vanishes quadratically at the Dirac point. Disorder is known to be perturbatively irrelevant, and previous theoretical work has assumed that the Dirac semimetal phase, characterized by a vanishing density of states, survives at weak disorder, with a finite disorder phase transition to a diffusive metal with a non-vanishing density of states. In this paper we show that nonperturbative effects from rare regions, which are missed by conventional disorder-averaged calculations, instead give rise to a nonzero density of states for any nonzero disorder. Thus, there is no Dirac semimetal phase at non-zero disorder. The results are established both by a heuristic scaling argument and via a systematic saddle point analysis. We also discuss transport near the Dirac point. At the Dirac point, we argue that transport is diffusive, and proceeds via hopping between rare resonances. As one moves in chemical potential away from the Dirac point, there are interesting intermediate-energy regimes where the rare regions produce scattering resonances that determine the DC conductivity. We derive a scaling theory of transport near disordered 3D Dirac points. We also discuss the interplay of disorder with attractive interactions at the Dirac point, and the resulting granular superconducting and Bose glass phases. Our results are relevant for all 3D systems with Dirac points, including Weyl semimetals.


I. INTRODUCTION
The discovery of two-dimensional (2D) Dirac systems such as graphene and the surface states of topological insulators has sparked an explosion of activity in condensed matter physics [1,2]. Such materials, which are gapped everywhere except at isolated points in the Brillouin zone, play host to an abundance of new physics. In particular, when the chemical potential is placed at the 'Dirac point', they display behavior that is intermediate between metals and insulators, in that the spectrum is gapless, but displays a vanishing low-energy density of states (DOS). The theoretical prediction [3][4][5][6][7] and experimental discovery [8][9][10][11] of three dimensional (3D) Dirac points provides a higher dimensional version of this behavior, and has ignited a blaze of interest in 3D Dirac points.
When the Fermi level lies precisely at the Dirac point of a clean 3D system, the density of states vanishes and the mean free path diverges. The consensus in the theory literature, from original work by Fradkin [12] in the 1980's to more recent work on Weyl semimetals [13][14][15][16][17][18][19], is that weak disorder is perturbatively irrelevant at 3D Dirac points, so that sufficiently weak disorder does not affect the vanishing of the density of states (DOS) at the Dirac point, or the divergence of the mean free path. Thus it has been assumed that there is a 'Dirac semimetal' phase (characterized by a vanishing DOS), which survives at weak disorder, and which undergoes a quantum phase transition at a critical disorder strength to a diffusive metal.
In this paper we show that this longstanding theoretical consensus is inaccurate, and that the Dirac semimetal phase does not exist at non-zero disorder. The source of the inaccuracy is non-perturbative rare region effects, which can dominate the physics at particle-hole symmet-ric points [20,21], and which were ignored in all previous analyses. When these rare region effects are correctly accounted for, the density of states at the Dirac point remains non-zero even for arbitrarily weak disorder, and the mean free path remains finite. There is no quantum phase transition at finite disorder. Rather, the Dirac semimetal only exists in the limit of vanishing disorder strength.
There is some similarity between the rare region effects discussed in this paper and the phenomenon of Lifshitz tails [22][23][24][25], of which a remarkably clear exposition can be found in [26]. However, there are also important differences. Whereas Lifshitz tails involve exponentially localized states [27] which exist inside a band gap, in the problem of interest to us there is no band gap, and thus a straightforward mapping to the Lifshitz tail problem is clearly impossible. The origin of the non-zero DOS in the present problem is more subtle, with the non-zero DOS arising due to power law bound resonances which co-exist with an extended continuum. The spinor nature of the wave function is also an essential requirement for the rare region effects we analyze, which are thus particular to Dirac fermion systems.
Our paper is organized as follows: in Section II we introduce the model of interest to us, and explain why previous theoretical analyses have (plausibly but erroneously) concluded that weak disorder can be ignored. In Section III we provide a heuristic scaling argument that suggests that a nonzero density of states arises for any nonzero disorder. In Section IV, we re-derive this result by means of a systematic saddle point analysis, conclusively establishing that the density of states is nonzero for any nonzero disorder. In Section V, we examine transport near the Dirac point, and discuss the multiple distinct transport regimes that arise as we tune the chemical potential away from the Dirac point. In Section VI we discuss the interplay of disorder and interactions, paying particular attention to the granular superconductor and Bose glass phases that may arise. We summarize our results in Section VII. The Appendix derives some results on massless three dimensional Dirac equations in spherically symmetric potentials. These results are used extensively in the main text.

II. MODEL AND BACKGROUND
The low energy Hamiltonian of interest takes the form where the two-component spinor ψ a (k) represents a state near the Dirac point a, with a momentum k relative to the Dirac point, ψ a (r) is its Fourier transform, and V (r) is a random scalar potential which is short range correlated and has mean zero. In general in a condensed matter system the dispersion about each Dirac point would be anisotropic, but for simplicity we consider the isotropic case.
In any lattice model with emergent Dirac fermions, Dirac points always come in pairs. For simplicity, we put a UV cutoff on the random potential so that it does not produce scattering between Dirac points, thus the different Dirac points are all decoupled. It is therefore sufficient for us to consider a single Dirac point with random scalar potential disorder, i.e.
This model most clearly exposes the relevant physics. Additional Dirac points can be retained in the analysis without changing the essential results. Near the Dirac point, in the clean limit V = 0, the low energy DOS (per Dirac point per unit volume) then vanishes as 3 . We now add weak quenched scalar potential disorder (strong disorder has been studied in [28]). A simple and highly intuitive argument for the irrelevance of scalar potential disorder proceeds as follows. An energy scale E sets a length scale v/E. Assuming short range correlated disorder with V (r) = 0 and V (R)V (R + r) = µ 2 0 f (r/b), with f (0) = 1 and f (x) decaying exponentially for x > 1, and averaging the disorder over a volume ( v/E) 3 using the central limit theorem, we conclude that the magnitude of the average potential over a length scale v/E will be |δV | ∼ µ 0 E 3/2 . The ratio |δV |/E ∼ µ 0 E 1/2 vanishes as E → 0, so one might conclude that in the asymptotic zero energy limit, the typical average potential vanishes more rapidly than the energy itself, and can thus be ignored.
An alternative argument for the perturbative irrelevance of disorder proceeds [14,29] by evaluating the electron self energy Σ, which yields Σ(ω, This vanishes more rapidly than ω at low energies and thus allows existence of sharp quasiparticles. Similarly, a self consistent Born approximation (SCBA) for the mean free path l leads to at the Dirac point; Λ is a UV cutoff and ν(E) ∼ E 2 . For sufficiently weak disorder µ 0 → 0, this admits only the trivial solution 1/l = 0 [at energies away from the Dirac point, l diverges as l ∼ 1/(µ 2 0 b 3 E 2 ) within SCBA]. This is in sharp contrast to two-dimensional Dirac materials, where, within SCBA, disorder produces a crossover to diffusive behavior at long length scales [30]. The difference arises because the DOS vanishes more rapidly in 3D, making disorder perturbatively irrelevant instead of marginal.
We note that in the above equation, we introduced a UV cutoff on the Dirac equation Λ. This in turn defines a length scale a = v/Λ, which should be of order the lattice scale. We henceforth set = v = a = 1 for convenience. All lengths are measured in units of a, all energies are measured in units of v/a and all times are measured in units of a/v. This consistent set of natural units will be used throughout this paper, although on occasion we will choose to display the factors of , v and a explicitly.
Finally, a simple RG analysis also suggests that scalar potential disorder is irrelevant.
The argument proceeds as follows: working with the Matsubara field integral and ensemble averaging over Guassiandistributed disorder using the replica trick (see e.g. [20]) gives rise to a quartic term of the form ∼ µ 2 0 dτ dτ d 3 xψ(x, τ )ψ(x, τ )ψ(x, τ )ψ(x, τ ). Straightforward power counting then reveals that µ 0 is irrelevant in the renormalization group sense at tree level.
As a result of all these excellent and intuitive arguments for the irrelevance of disorder, it has long been believed that there exists a Dirac semimetal phase characterized by a vanishing density of states and a diverging mean free path, which survives at weak disorder, with a quantum phase transition to a diffusive metal occurring at a nonzero critical disorder strength. It is the objective of this paper to establish that this belief is incorrectthere is no Dirac semimetal phase at non-zero disorder.
Rather the system has a non-vanishing density of states and a finite mean free path for any nonzero disorder. This new result arises due to the effect of exponentially rare regions which host critically localized resonances, and which are missed by conventional disorder-averaged calculations. An analysis of these rare regions effects will occupy the majority of this paper.

III. A HEURISTIC ARGUMENT FOR THE IMPORTANCE OF DISORDER
In this section, we provide a heuristic argument for the importance of disorder. To this end, we will invoke a somewhat artificial model of disorder, which nevertheless captures the essential physics (a more realistic treatment of disorder will be provided in the following section). The model of disorder we consider is one where we introduce Poisson distributed 'impurities', with a mean density of impurities n. Each impurity consists of a spherically symmetric scalar potential of the form , where b is a fixed length scale, λ (the 'strength' of the impurity) is taken from a Gaussian distribution P (λ) with mean zero and variance µ 2 0 /(nb 3 ), and ε obeys ε(b) = 1 and falls off at least as fast as 1/r 4 at long distances. i.e. each impurity is modeled as a scalar potential well with a 'tail' that falls off at least as fast as 1/r 4 . The precise form of the tail will not be important for our argument. We work in the weak disorder limit, which corresponds to µ 0 1/b and µ 2 0 nb. For convenience, we now take nb 3 = 1, although our results are readily generalized to other cases with no important changes.
It can be shown (see Appendix) that a single impurity will trap a bound state for specific values of λb. The bound states take the form where f and g are purely radial functions with no angular dependence, and the φ ± are two component spinors (detailed in the Appendix) which have definite total angular momentum j and which have orbital angular momentum differing by one. Bound states arise with all values of j. However, bound states with high j require a deeper or wider well than bound states with smaller j. In the weak disorder limit, the physics near the Dirac point is dominated by bound states on rare very strong impurities with j = 1/2 (which are linear superpositions of states with orbital angular momentum zero and one). We emphasize that for the j = 1/2 solutions, |ψ| 2 is isotropic i.e. the j = 1/2 bound state solutions have an isotropic probability density. The existence of bound states is not sensitive to the form of ε(r), and it is convenient to take an ε(r) that falls off infinitely fast (corresponding to a square well). In this case the bound states arise at λ = λ c ≈ mπ/b, for nonzero integer m, and have f (r → ∞) ∼ 1/r 2 and g(r → ∞) = 0. The position of λ c shifts if we make a different choice for ε(r), and the subleading piece g(r) changes (e.g. g(r) ∼ 1/r 5 for ε(r) ∼ 1/r 4 ) but the existence of bound states and the leading 1/r 2 falloff of the wave function do not change, and the scaling argument that we will now present is unaltered. We therefore stick to the 'square well' potential, which has bound states for λb ≈ mπ. Recall that P (λ) ∼ exp(−λ 2 /2µ 2 0 ). In the weak disorder limit the physics is dominated by the j = 1/2 bound states, which have m = ±1 and proba-bility density P (λ c ) ∼ exp(−Cπ 2 /(2µ 2 0 b 2 )), whereC is a numerical prefactor of order one. Now, the bound states in question are rather delicate. Unlike the case of Lifshitz tails, where the bound states are exponentially bound, here the states are only power law bound. Moreover (also unlike the case of Lifshitz tails), true bound states arise only for precise values of λ = λ c . In an infinite sample, even an infinitesimal deviation from λ = λ c or E = 0 leads to the disappearance of the bound state (indeed the solution becomes nonnormalizable). However, the disappearance of the bound state becomes apparent only at a very large length scale R, which diverges as λ → λ c and E → 0. On length scales less than R, the solution for λ close to but not quite λ c is indistinguishable from the truly bound state. In a system with a nonzero density of impurities hosting almost bound states, where R exceeds the spacing between such almost bound states, such small deviations from λ c should not matter. Thus, we expect that the density of quasi-bound states will be where δλ is a to-be-determined quantity that tells us how close we have to get to λ c in order to have a state that looks effectively bound. Moreover, these quasi-bound states will not all be strictly at zero energy. Rather they will be spread over an energy window of width δE. The contribution to the density of states coming from quasi bound states will then take the form where we have defined To make further progress requires understanding what happens when λ and E are slightly perturbed from λ c and 0 respectively. To this end, it is instructive to calculate the scattering cross section σ(λ, E) of a single well. The scattering cross section σ is given by the formula where δ is the phase shift [44]. We now have to determine the phase shift δ. Since we have already established that the nature of the 'tail' of ε(r) does not qualitatively alter the physics, we model the well as being simply a square well -this greatly simplifies the calculation.
Application of standard trigonometric identities, as well as the results from the Appendix for A and B (see also [45]), then leads to the result This equation contains a great deal of physics. The magical values of λ = λ c which give rise to bound states are revealed as resonances, which correspond to phase shifts δ = π/2. These resonances 'pull' some density of states out of the continuum and down to zero energy (in the form of bound states). Meanwhile, in the scaling limit E → 0, the cross section is a tightly peaked Lorentzian, with From this we conclude that there is a line of resonances in the λ, E plane, and that these resonances have width ∼ bE 2 in both λ and E. This leads us to the scaling δλ ∼ δE ∼ bE 2 . Substituting this into (6) tells us that the low energy density of states is just ν 0 , given by (7). When making this estimate, we have not taken the non-resonant 'extended' states into account, thus this estimate is valid only on scales E < √ ν 0 , where the density of states ν 0 from special wells exceeds the DOS ∼ E 2 from the extended states. However, in this regime, our scaling theory reveals that the DOS is given by (7), and is non-zero for arbitrarily weak disorder. This establishes that even though disorder is perturbatively irrelevant at the Dirac point, it cannot be neglected. We close this section by highlighting one important point. The white noise limit b → 0 is 'pathological'. In the limit b → 0 at any, even very small, fixed µ 2 0 b 3 (necessary to keep the disorder strength constant), the system can not remain in the fully weak disorder limit µ 0 1/b and these 'rare quasi bound states' become not at all rare. This aspect of the white noise limit will become clear in the saddle point calculation presented in the following section.

IV. A SYSTEMATIC CALCULATION OF THE DENSITY OF STATES
In the previous section we provided a heuristic scaling argument suggesting that the DOS is non-zero in the presence of arbitrarily weak disorder. We now rederive this result using 'standard' techniques. We follow the route taken in [26], suitably generalized to the present problem. We begin by noting that we are dealing with a system governed by the Dirac Hamiltonian where σ i is a Pauli matrix, and the ψ V n are the twocomponent spinor eigenfunctions which satisfy the above equation for eigenenergies E V n in a given random scalar potential V . Repeated indices are summed over. The density of states per unit volume at an energy E and for a given disorder configuration V , ν V (E), can be expressed as where L is the linear size of the system. We now introduce a spinor Lagrange multiplier field χ and a scalar Lagrange multiplier Υ to rewrite this as Integrating out the Lagrange multiplier fields gives us a delta function which picks out only those configurations ψ(x) which satisfy the Dirac equation with E = E V n and are properly normalized. The functional integral over ψ(x) then reproduces (12). Note that in order for the exponent to be properly dimensionless when ψ is properly normalized, the scalar Υ must be a pure number, whereas So far we have discussed the DOS for a specific disorder realization. We now average over disorder (assuming that the disorder is a Gaussian random variable, which is short range correlated with a correlation length ξ), to obtain a disorder averaged density of states ν, which takes the where We have performed one formal manipulation, defining the rescaled variablesχ † = iχ † W andΥ = iΥ. We have scaled the lengths by the microscopic length scale a (which has been set equal to one), soχ † and ψ are now dimensionless. W measures the disorder strength, with W 2 ∼ µ 2 0 ξ 3 . Meanwhile, K is the correlation function for the disorder, and we have defined K −1 according to . We assume that K is an isotropic and normalized function which is short ranged with a characteristic scale ξ (for definiteness, we could take K to be a normalized isotropic Gaussian with width ξ, but the results will be independent of the precise shape of K).
We now make the one essential approximation required by our approach: we calculate the density of states ν in a saddle point approximation. The saddle point equations obtained by varying V ,χ,Υ, ψ † , and ψ respectively are: Now, we note that the equation forχ † , (20), is just the Hermitian conjugate of the equation for ψ, (17). Thus, we takeχ † = χ 0 ψ † , where χ 0 is a scalar. We want to search for solutions at real energies, so we want the saddle point Hamiltonian to be Hermitian. This then demands that we should take χ 0 to be a real number, although it could be either positive or negative. Substituting into (16) tells us that within the saddle point approximation, the DOS is dominated by potential configurations with Thus we find that the equations of motion all boil down to a single (non-linear) integro-differential equation which takes the form (21) where we have specialized to E = 0 and require that the solutions ψ be properly normalized d 3 x|ψ(x)| 2 = 1. We emphasize that we are allowed to tune χ 0 in order to find a solution. The contribution of a particular solution to the disorder averaged density of states is found by substituting the saddle point solution into (14). This yields We must sum over all saddle point solutions to accurately obtain the density of states. We note that for any solution ψ 0 of the above equation with E, χ 0 , there will be a corresponding solution to the same equation with E → −E, χ 0 → −χ 0 , ψ 0 → Cψ 0 , where C is the particle hole symmetry operator, and this corresponding solution will have the same cost action. Thus, the DOS will be even in energy. It is sufficient for our present purposes to determine the DOS at E = 0.

A. SCBA
The equation (21) clearly has saddle point solutions that are plane waves ∼ 1 L 3/2 e ik·x , with k = |E|. Substituting into (22) tells us that the contribution of the plane wave saddle points to the density of states is Thus, the density of states is just equal to the number of plane wave solutions at a given energy, divided by the volume. In the E → 0 limit, the number of plane wave solutions in a window of energies between E and E + δE scales as L 3 E 2 δE, thus, the density of states coming from these 'plane wave' saddle points scales as E 2 , and vanishes at zero energy. The various arguments for the perturbative irrelevance of disorder outlined in Section II essentially amount to the statement that perturbation theory about these translation invariant saddle points converges. However, we will now proceed to show that there are additional solutions which satisfy the saddle point equations, and while these other solutions have a weight that is exponentially small in weak disorder, the density of these other saddle points does not vanish as E → 0. It is these other saddle points (which correspond to rare resonances) this will give rise to a non-zero density of states at E = 0.

B. A saddle point treatment of rare regions
Motivated by the scaling analysis in the previous section, we look for localized and normalizable solutions to the saddle point equations. We assume for convenience that the localized solution is centered at the origin. It is convenient to interpret (21) as a Dirac equation in an effective potential Guided by our previous analysis, we expect to find a solution of the form (4) with total orbital angular momentum j = 1/2, and such a solution has the property that |ψ| 2 is isotropic. If K is an isotropic function, it then follows that the effective potential is spherically symmetric. Thus, (21) can be re-interpreted as a Dirac equation in a spherically symmetric potential, where the shape of the potential must be determined self consistently, and where the potential profile tracks the probability density for a bound state wave function with total angular momentum j = 1/2. (There will also be solutions corresponding to higher values of j, but these will require a deeper or wider self consistent potential well and thus will have exponentially suppressed contribution to the DOS, such that the DOS will be dominated by solutions with j = 1/2). Now, we have already assumed that the kernel K is sharply peaked e.g. K(r) ∼ 1 ξ 3 exp(−r 2 /ξ 2 ). On long length scales r ξ, K can be modeled as a delta function, leading to the simpler equation Making the ansatz (4) with j = 1/2, this can be written out in components as where we have introduced the shorthand variable X = χ 0 W/4π. We now imagine constructing a power series where each successive term is higher order in 1/r. Guided by our earlier work on the linear Dirac equation (see Appendix), we look for a solution where f ∼ 1/r 2 and g falls off faster. This leads to a long distance solution that has the form with an undetermined scale factor A. We can readily see that this takes the form of a perturbation series in the small parameter X 2 A 4 /r 6 1. We now check for self consistency. We note that the solution identified above describes an effective potential that falls off as 1/r 4 at large r. We have already identified (see Appendix) that the Dirac equation in a spherically symmetric potential that falls off as 1/r 4 has a solution where f ∼ 1/r 2 and g ∼ 1/r 5 . Thus, the solution we have constructed is indeed a correct self consistent solution of the non-linear integro-differential equation (21) at large distances.
The expansion introduced above breaks down at r c = X 1/3 A 2/3 = (χ 0 W A 2 /4π) 1/3 . A numerical investigation of the equation (25) reveals that solutions with the asymptotics identified above are generally singular at r = 0. However, at small distances r ≤ ξ, modeling the disorder correlation function K as a delta function is clearly inappropriate, and thus we cannot work with (25), but instead we must work with the full integrodifferential equation (21). On distances r < ξ, the convolution with K produces an effective potential which is roughly constant.
This then implies that we are solving the Dirac equation in a spherically symmetric well that is of constant depth ∼ χ 0 W on length scales less than ξ, but has a 1/r 4 tail. We know what the solutions to this problem look like (from the Appendix): they have |ψ| ≈ A at short distances r < ξ and |ψ| ≈ Aξ 2 /r 2 at long distances r > ξ. Normalization fixes A 2 ≈ ξ −3 . The self consistent potential defined by (24) then is uniform with depth χ 0 W A 2 = χ 0 W ξ −3 at short distances r < ξ and falls off as χ 0 W ξ/r 4 at long distances r > ξ. We recall that we are free to tune χ 0 to find a solution. Given the results derived in the Appendix, we expect that bound states will exist for an infinite discrete set of χ 0 . However, since the contribution to the DOS falls off exponentially with χ 0 (Eq. 22), the dominant contribution will come from the smallest value of χ 0 that allows us to have a solution. The smallest value of χ 0 that allows for a solution has χ 0 W A 2 ξ ≈ 4 (see Appendix, Fig. 2) i.e. χ 0 ∼ 1/(W ξA 2 ).
Thus we have shown that there exists a normalizable and localized solution to the saddle point equations where the wave function falls off as A/r 2 at large distances, with A 2 ∼ ξ −3 and χ 0 ∼ 1 W ξA 2 . Substituting this into the expression for the DOS (22) and approximating we find that the saddle point solution identified above makes a contribution to the DOS of order (29) where in the final expression we have restored , v and a for clarity. Here C is a numerical pre factor (expected to be of order one), which cannot be determined without actually solving the full integro-differential equation (as opposed to showing a solution exists and identifying its asymptotics).
We note that thus far we have evaluated the contribution to the disorder averaged DOS from a single localized saddle point solution centered at the origin. However, the localized saddle point solution could be centered anywhere in the sample and we must sum over all possible locations of the bound state centre i.e. there are an extensive number of saddle points of this form contributing to the density of states. The summation over the centre of mass co-ordinates cancels the 1/L 3 factor in the above equation. Another way to state this result is to note that while the density of plane wave solutions to the saddle point equations vanishes as E → 0, the density of localized solutions to the saddle point equations does not vanish as E → 0, so that the localized solutions to the saddle point equations actually dominate the low energy DOS. This we obtain a final expression for the contribution to the DOS from localized solutions to the saddle point equation that takes the form Determining the precise constant C in the exponential requires determining the precise shape of the localized solution everywhere (i.e. not just the asymptotics), whereas determining the pre-exponential factor requires a consideration of fluctuations about the saddle point (for more details on this procedure, see [26]). We defer consideration of these issues to future work. However, we note that the scaling of the DOS with disorder strength and well radius (exp(−ξ/W 2 ) ) is the same as that from the heuristic scaling approach employed in Section III, if we identify W with the rms potential in the well W 2 a 3 ∼ µ 2 0 b 3 , and if we identify b ∼ ξ. Thus, a systematic saddle point calculation reveals that the heuristic scaling approach developed in Section III obtains essentially correct results. It also reveals the flaw in the perturbative arguments detailed in Section II: those arguments only consider fluctuations about the wrong (i.e. translation invariant) saddle point, whereas the physics is dominated by different, translation symmetry breaking, localized saddle points.
We note that we have only taken into account the localized solutions to (21) which minimize the cost action (the solutions with j = 1/2). There will be additional localized solutions with higher total angular momentum, but these will have a larger cost action, and hence will make an exponentially smaller contribution to the DOS. Still, given the likely existence of higher angular momentum saddle points, (30) should properly be viewed as a lower bound on the DOS.
We note that in the white noise limit ξ → 0, the contribution from these localized solutions become of order one. However, the white noise limit is pathological, for the following reason. The saddle point solutions that give rise to density of states at E = 0 involve potential fluctuations of magnitude W χ 0 A 2 ≈ W 1 W ξA 2 A 2 ≈ 1/ξ. In the white noise limit the bound states require increasingly large potential fluctuations. If we work with a model of unbounded disorder (such as the model used in this section) then we obtain an order one density of states in the limit ξ → 0. However, in this limit the solutions are singular at r → 0, and require potential fluctuations of diverging magnitude. If the disorder fluctuations are ultimately bounded, then a different analysis is required. For disorder fluctuations that are bounded by cΛ, the analysis in the present subsection applies for ξ > a/c. Meanwhile, the weak disorder limit is W 2 /ξ 1, or equivalently µ 0 b 1.

V. TRANSPORT NEAR THE DIRTY DIRAC POINT
Thus we have shown that (notwithstanding perturbative arguments to the contrary), the density of states at a disordered 3D Dirac point does not vanish, even for arbitrarily weak disorder. We now turn our attention to the transport properties. A systematic approach to transport properties would involve writing down a supersymmetric sigma model (or a replica sigma model), and incorporating the effect of the localized saddle point solutions identified in Section IV. (For a discussion of how to translate a saddle point calculation of the form developed in Section IV to the supersymmetric and replica formalisms, see [26]). However, emboldened by the success of our scaling arguments in calculating the density of states, we now choose instead the simpler and more intuitive option of generalizing our scaling arguments to a scaling theory of transport near the 3D Dirac point.
That will be the focus of this section.
We note that while traditional Lifshitz tails involve exponentially bound states that live in a band gap, here we are dealing with power law bound resonances that co-exist with a continuum of extended states (albeit a continuum that has vanishing density of states). The resulting transport behavior will be very different to that encountered with traditional Lifshitz tails.
We recall that a single rare potential well with width b and depth λ has a cross section for states at an energy E that takes the form i.e. there is a line of resonances in the (λ, E) plane, with width δλ ∼ δE ∼ bE 2 . It is instructive to calculate the mean free path from scattering off resonant rare regions. This behaves as where ν 0 is given by (7) and we recall that we are working with a model of disorder where P (λ) ∼ exp(−λ 2 /2µ 2 0 ). At high energy where the SCBA remains valid, the resulting mean free path is l ∼ 1/(µ 2 0 b 3 E 2 ). The rare regions start to dominate the scattering when this SCBA mean free path exceeds that due to the rare regions, which is at an energy scale E < ∼ √ ν 0 /(µ 0 b). However, rare regions do not start to dominate the density of states until E < ∼ ν 1/2 0 (which is a much smaller energy scale, in the weak disorder limit µ 0 b 1). Moreover, we do not enter the strong scattering / hopping conduction regime until E < ν 0 b (according to the Ioffe-Regel criterion [46]). Thus, we are led to identify four distinct regimes. At the highest energies E ≥ √ ν 0 /µ 0 b, the behavior is governed by SCBA. For ν 1/2 0 < E < ν 1/2 0 /µ 0 b, the DOS is dominated by extended states, but the (still weak) scattering is dominated by rare regions. Meanwhile, in the regime ν 0 b < E < ν 1/2 0 , the DOS and scattering are dominated by the rare regions, but the mean free path is still much longer than the wavelength and the scattering is in this sense weak. Finally, for E < ν 0 b, the mean free path is less than 1/E, and we are in the 'strong scattering' regime where it no longer makes sense to talk about weakly-scattered extended states. In this regime, the states all live on rare regions, and transport proceeds by hopping. In this region we have δλ ∼ δE ∼ bE 2 ∼ b 3 ν 2 0 , and the typical hopping rate is also b 3 ν 2 0 . Meanwhile, the density of rare regions is P (λ c )δλ ∼ b 3 ν 3 0 , and the typical spacing is (bν 0 ) −1 . Thus, transport in this regime occurs due to hopping over length scales (bν 0 ) −1 .
In both intermediate energy regimes, the carriers spend a typical time ∼ b −1 E −2 trapped on each resonant special well (this is just the width of the resonance) whereas the time spent traveling freely in between special wells is proportional to the mean free path l ∼ (bν 0 ) −1 . Thus, in the intermediate energy regime ν 0 b < E < ν  Table I. In each case the diffusivity is D ∼ l 2 /τ , with l the typical hopping distance in the hopping regime and the mean free path in the other regimes. The time between hops or scattering events is τ . The zerotemperature conductivity for these noninteracting carriers is then σ DC = νe 2 D, where ν is the DOS. Stitching together the low energy (hopping dominated) and high energy (SCBA) regions leads to the plot Fig.1.
We note that when estimating the diffusion constant we ignore the possibility of interference between distinct paths. Such interference could give rise to localization or anti localization behavior in the hopping model at very long length scales. Now the fact that the wave functions have a 1/r 2 falloff (which is slower than 1/r d ) guarantees that there can be no localization on an individual resonance. However, at the very longest length scales we have a theory of non-interacting fermions hopping on a random network of (exponentially widely spaced) resonances, and we may worry about interference between distinct paths on this network. Within the model of purely scalar potential disorder considered here, the various Dirac points are all decoupled. It is widely believed that one cannot localize a single Dirac fermion. This belief is based on calculations involving sigma models, which are generally designed to treat fluctuations about a translation invariant saddle point. We have shown that a translation invariant saddle point is not the appropriate starting point, and thus the sigma models need to be rederived. Assuming that it remains impossible to localize a single Dirac fermion even taking rare region effects into account, the only possibility to be wary of is anti localization. However, the standard scaling arguments suggest that anti localization should be a weak effect in three dimensions (at least for weak disorder W → 0), with the β function for the conductance taking the form β(g) ∼ 1 + f (W ), where f (W → 0) = 0. Thus we conclude that the neglect of interference between distinct paths is not a real problem, and that transport at the lowest energies should indeed be diffusive, and dominated by hopping between rare resonances.

VI. INTERPLAY OF DISORDER AND INTERACTIONS
Thus far we have concentrated on non-interacting 3D Dirac fermions. We now turn our attention to the interplay of disorder and interactions. Above a critical interaction strength, repulsive interactions destroy the Weyl semimetal phase [31,32,47]. Subcritical repulsive interactions suppress (charged) rare regions, and reduce the rare-region DOS at the Dirac point. We defer further

Energy regime
Description Length scale Time scale DOS Diffusivity DC conductivity  Table listing  The rare regions dominate the scattering for all regimes other than the highest-energy SCBA regime. The nonzero slopes on this log-log plot are ±2. For more details, see text and Table I .
consideration of repulsive interactions to future work. Instead, we now consider the interesting interplay that occurs between disorder and attractive interactions at the Dirac point. Attractive interactions above a critical strength will trigger superconductivity in the clean system [33,34]. Subcritical interactions will produce local pairing on rare regions where the local DOS is non-zero over a larger length scale than the local coherence length Ξ. Establishment of phase coherence between islands by Josephson coupling will then drive the system into a (granular) superconducting state at sufficiently low temperatures. We have discussed similar phenomena for the 2D Dirac system in [20]. We focus on estimating the energy scale for the superconducting state, in the presence of a random scalar potential that is approximately Gaussian distributed (for small fluctuations), but which is ultimately bounded, with no local fluctuations that are larger than Λ.
Local pairing occurs in islands of local average potential µ and size L ≥ Ξ, where Ξ ∼ (v/ω D ) exp(1/Gµ 2 ) is the local coherence length in the BCS approximation, ω D is the Debye frequency and G is the strength of the attraction in the leading pairing channel. Integrating over L in a saddle point approximation, we find the result is dominated by islands of size L ∼ = Ξ. The probability of finding such an island is where G −1/2 marks the boundary of the weak coupling BCS regime. This is dominated by the regions near the upper limit of this integral, and yields Here G 1 = 1 Λ 2 , and G c is the critical coupling for superconductivity in the clean system. This density of superconducting islands is doubly exponentially small in G for G → 0 when even the maximally doped islands with local µ ≈ Λ have to be exponentially large; but is only exponentially small in G for intermediate G, when small superconducting islands with local doping µ ≤ Λ can form.
In the intermediate range of G, the energy scale for local Cooper pairing in each of the dominant islands is of order ω D . However, the sample will exhibit global superconductivity only if phase coherence is established between islands. The Josephson coupling between distant islands J may be determined by generalizing the calculation in [35] to the 3D Dirac point. We find that J ∼ 1/r 5 . Since the Josephson coupling falls off with distance faster than 1/r 3 , the coupling between nearest neighbor islands dominates. The system of locally superconducting islands embedded in a semimetal then establishes global phase coherence on temperature scales smaller than the typical nearest neighbor Josephson coupling. This leads to an estimated critical temperature for phase ordering We have implicitly assumed that the pairing is s-wave. If the 'local pairing' was not s-wave, then the Josephson couplings would be frustrated, and the ground state would be a 'gauge glass' [36]. We leave further discussion of non-s-wave orders to future work, noting only that in [33] it was determined that δ-function attraction in the clean system favors s-wave pairing.

A.
Attractive and repulsive interactions: We now discuss the situation when Coulomb repulsion coexists with retarded attractive interactions. We assume that the Morel-Anderson condition [37] is satisfied, so that local pairing on islands still occurs. However, the effective Hamiltonian for the islands must now contain not only the Josephson couplings, but also charging effects (electrostatic interactions may be neglected due to screening [38]). Thus, the effective Hamiltonian for the islands is where i and j label superconducting islands, φ i is the phase of the i th island, and n i = i∂/∂φ i . The Josephson couplings J ij operate primarily between nearest-neighbor islands, as previously discussed, and the V i n i term reflects the random scalar potential on the islands. Such Hamiltonians have been long discussed in the theory literature [39][40][41], and are known to support a superconducting phase, and also a Bose glass [42]. The glassy phase is characterised by an infinite superconducting susceptibility, but no long range order, and has a regime of stability that grows larger as the system becomes more disordered.

VII. CONCLUSIONS
Thus we have demonstrated that a 3D Dirac point has a non-vanishing density of states ν 0 ∼ exp(−ξ/W 2 ) for weak scalar potential disorder with strength W and correlation length ξ. The physics at low energies is dominated by exponentially rare, power law bound resonances which break translation symmetry and 'pull' density of states down to zero energy. We have shown how the density of states can be estimated using a rare regions scaling argument, and also using a systematic saddle point analysis. The systematic saddle point analysis also reveals what was missed by the existing theoretical arguments for the irrelevance of disorder (detailed in Section II): those arguments only considered fluctuations about a translation invariant saddle point, whereas the non-zero density of states arises due to other, translation non-invariant saddle points, which cannot be accessed through perturbation theory about a translation invariant saddle point.
We have also constructed a scaling theory of transport near a 3D Dirac point in the presence of random scalar potential disorder with strength W and correlation length ξ. This theory reveals that there are four distinct transport regimes. At the highest energies |E| > √ ν 0 ξ/W , the SCBA solution applies and both scattering and the DOS are dominated by extended states. The DOS scales as E 2 and the DC conductivity is constant. For √ ν 0 < |E| < √ ν 0 ξ/W , scattering is dominated by rare regions, but the DOS is still dominated by extended states. In this regime the DOS and the DC conductivity both scale as E 2 . For ν 0 ξ < |E| < √ ν 0 both scattering and the DOS are dominated by rare regions. The DOS is constant, but the conductivity scales as E 2 . Finally, for |E| < ν 0 ξ, we enter a 'strong scattering' regime in which we argue that both the DOS and DC conductivity saturate to non-zero constants. Finally, we have also discussed the interplay of attractive interactions with rare resonances, which can drive the system into a granular superconducting phase, with a critical temperature that we estimate. We have also discussed the Bose glass phases that can arise in the presence of both attractive and repulsive interactions.
This work has established that the existing framework for thinking about 3D Dirac points -in terms of translation invariant disorder averaged theories -is inaccurate at the lowest energies. Instead, one must take into account the effects of rare resonances, which control the physics close to the Dirac point. In light of the rapid experimental advances in synthesizing materials supporting 3D Dirac points, we hope that it will soon be possible to probe the asymptotic low energy regime in experiments, and to directly test the scaling theory advanced in this paper.

VIII. APPENDIX
In this appendix we solve the three dimensional Dirac equation in a spherically symmetric potential. The results obtained in this way will be essential to construction of our argument. The Dirac equation in a spherically symmetric potential can be written as We work with a single Dirac point, since the different Dirac points are all decoupled. The resulting equation for a two component spinor wave function is sometimes also referred to as the Weyl equation. However, we continue to refer to it as a Dirac equation here, to emphasize that our results are not particular to Weyl semimetals. The eigenstates of the Dirac Hamiltonian are also eigenstates of total angular momentum j, but are not eigenstates of orbital angular momentum l. Using the standard Pauli matrix multiplication identity σ i σ j = δ ij + iε ijk σ k , we rewrite the gradient term as using the notationr = r/r and r 2 = r j r j , and where L is the usual quantum mechanical angular momentum operator. This prompts us to search for a solution of the form ψ = R(r)φ, where R is a scalar function that depends purely on radius, whereas φ is a two component spinor which is an eigenstate of the angular momentum operator, and which is independent of radius. Now, the eigenstates of the operator σ · L are two component spinors φ ± j,jz with total angular momentum j, angular momentum projection onto the z-axis j z , and orbital angular momentum l ± = j ∓ 1/2, which take the explicit form [43] where the Y functions are the usual spherical harmonics. We note that the ± superscript refers to the angular structure. Using the identities J = L + 1 2 σ and J · J = j(j + 1) 2 , L · L = l(l + 1) 2 , we can show that the spinors obey σ · Lφ ± j,jz = −(1 + κ) φ ± j,jz , where κ = −(j+1/2) is a negative integer for φ + and κ = j+1/2 is a positive integer for φ − .
We note that the functions φ ± j,jz have orbital angular momentum differing by one, and thus have opposite parity under inversion. Since σ ·r commutes with the angular momentum operator and changes sign under inversion, it follows that it must turn φ + into φ − and vice versa. Since the gradient term mixes the angular sectors φ ± , the eigenstates of the Hamiltonian must be linear superpositions of pieces with φ + and φ − angular structure. Thus we find that the eigenstates in the vicinity of the Dirac point take the form where f and g are purely radial functions with no angular dependence. Substituting this expression for the wavefunctions into the Dirac equation leads to the two equations where κ is a positive integer for one solution, and κ is a negative integer for its degenerate partner which differs only in its angular structure. Let us pick positive κ for specificity.

A. Square wells
We begin by considering a 'square well' potential V (r) = λΘ(b − r), although we will relax this approximation in due course. We note that because of the particlehole symmetry of the problem positive and negative V must yield identical results. Substituting the square well potential into (40) and performing some elementary manipulations then leads to the equation We recognize this as the spherical Bessel equation, whose solutions are spherical Bessel functions. Substituting f back into the equation for g then determines g. Thus, the solutions for arbitrary E = V take the form where J and K are Bessel functions of the first and second kind respectively. To save writing, we now adopt a system of units where v = 1. We will re-introduce v whenever necessary for clarity.
For r < b, V = λ. In this region, we must have B = 0 to have a regular solution at the origin. Meanwhile, for r > b, V = 0. In this region we can have A = 0 and B = 0. Thus, we have Since we are dealing with a first order differential equation, only the wave function need be continuous (there is no requirement that derivatives be continuous). Imposing continuity of the wave function then implies that where ∆ is the determinant of the 2 × 2 matrix. This fails for special values of E where the matrix is singular (vanishing determinant). We note that continuity of the wave function also implies continuity of the probability density (given by the norm squared of the wave function). The norm squared of the wave function at r = b (defined as |f | 2 +|g| 2 ) never vanishes, and scales as (λ − E) −2 b −2 in the limit of large |λ − E|b while saturating to a constant in the limit of small |λ − E|b. Thus, the probability density just outside the well never vanishes, and there is always 'leakage' of the probability density out of the region r < b. Moreover, the spherical Bessel functions only decay as 1/r at long distances, so the probability density only decays as 1/r 2 at long distances. Thus, the solutions constructed above are not normalizable in an infinite volume.
A qualitatively different (and properly normalizable) exterior solution exists when E = 0. When E = 0 then the two equations in (40) decouple for r > b, and can be straightforwardly solved to give an exterior solution Recall that κ is a positive integer. This corresponds to a bound state if and only if we pick the solution g(r) = 0, which comes about if g(r) is matched to a node of the interior Bessel function. This in turn happens only for special values of the well depth λ c . The probability density in this bound state decays like 1/r 4 outside the well (i.e. most of the probability density is localized on the well and the solution is properly normalizable). Although there is a well depth corresponding to a bound state for all values of κ, larger values of κ require a deeper (or wider) well in order to have a bound state. The physics of interest to us will thus be controlled by the minimal well, which has a bound state for κ = 1. Bound states with κ = 1 arise when λb ≈ mπ, where m is a positive integer. Again, values of m greater than one involve deeper or wider wells, and the minimal well which controls the physics has κ = 1 and m = 1, with a well depth λ c ≈ π/b. Note that there is a single parameter that must be tuned to get a bound state: either we can fix b and tune λ, or we can fix λ and tune b.
We note that the angular eigenfunction φ − has total angular momentum j = 1/2 (for κ = 1), but may have j z = ±j. Thus there are two bound states corresponding to the κ = 1 solution identified above. We note that there are two additional bound states corresponding to κ = −1 and λ = λ c = π/b, which now corresponds to an angular eigenfunction φ + and has f (r) = 0. Thus, there are four bound states per Dirac point for each special well.

B. Beyond square wells
Thus far we considered square well potentials. Now we consider a potential that has a long range tail. For specificity, we consider the potential V (r) = λΘ(b − r) + ε(r)Θ(r − b), where ε(r) = λb 4 /r 4 . The equations (40) for zero energy states in the domain r > b then become λb 4 r 4 f (r) = ∂ r g(r); − λb 4 r 4 g(r) = ∂ r + 2 r f (r).
It can be readily checked by solving the radial equations numerically on Mathematica that the 1/r 4 potential is not special. Bound states arise also for exponential tails, and for power law tails where the potential falls off faster than 1/r 4 . In all cases obtaining a properly continuous and normalizable bound state solution requires tuning a single parameter λb. One way to see that there is a single parameter which has to be tuned is the following: the interior solutions have an overall scale factor A. The exterior solutions have the overall scale factor C. Matching g fixes the ratio of scale factors A/C, but we still have to match f . Matching f requires tuning one parameter, and the relevant parameter here is λ 0 b. Moreover, in all cases the probability density decays as 1/r 4 at large distances, just as for the square well i.e. the asymptotic behavior is unchanged.
Thus, we have demonstrated that the square well potential is not special and that qualitatively similar be-