Breakdown of universality in three-dimensional Dirac semimetals with random impurities

Dirac-Weyl semimetals are unique three-dimensional (3D) phases of matter with gapless electrons and novel electrodynamic properties believed to be robust against weak perturbations. Here, we unveil the crucial inﬂuence of the disorder statistics and impurity diversity in the stability of incompressible electrons in 3D semimetals. Focusing on the critical role played by rare impurity conﬁgurations, we show that the abundance of low-energy resonances in the presence of diluted random potential wells endows rare localized zero-energy modes with statistical signiﬁcance, thus lifting the nodal density of states. The strong nonperturbative effect here reported converts the 3D Dirac-Weyl semimetal into a compressible metal even at the lowest impurity densities. Our analytical results are validated by high-resolution real-space simulations in record-large 3D lattices with up to 536 000 000 orbitals.


I. INTRODUCTION
The discovery of Dirac and Weyl semimetals (DWSMs) has provided a rich arena for probing novel gapless phases of matter with unique transport properties and topological features [1]. Several types of gapless systems featuring Dirac or Weyl points in three-dimensional (3D) momentum space have been realized [2][3][4]. The simplest DWSMs exhibit twofold or fourfold degenerate linear-band touching points at the Fermi level with isotropic velocities and a possible replication into disjoint momentum-space valleys. Their pointlike Fermi surface is protected against band gap opening due to either topological constraints-in Weyl systems with broken time-reversal (T )o ri n v e r s i o ns y m m e t r i e s( P) [ 1]-or crystal symmetries in TP-symmetric Dirac systems [5,6]. Thus any clean DWSM is an incompressible electron gas with a quadratically vanishing density of states (DoS). Despite the inefficient charge screening at the node, this paradigm is believed to survive electron-electron Coulomb interactions, giving way to a marginal Fermi liquid behavior [7,8]. In addition, weakly interacting DWSMs can display strongly renormalized Fermi velocities, but the nodes' integrity and topology are expected to remain robust [9][10][11][12][13][14][15][16].
An outstanding question is whether random on-site potentials ubiquitous in realistic systems (e.g., due to impurities in the crystal lattice) can give way to a compressible diffusive metallic phase with a finite nodal DoS [6,[17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]]. An early result by Fradkin predicted that Dirac nodes are stable in d = 2 + ǫ dimensions, below some critical disorder strength [17]. The robustness of DWSMs against weak random perturbations is best visualized by considering a massless particle moving through a short-ranged random potential of strength W . Since, near a node, the de Broglie wavelength, λ =hv/E , largely exceeds the disorder correlation length, the central limit theorem applies, and the fluctuations around the average potential inside a volume λ d must scale as δV ∝ W λ −d/2 .I n d = 3, the fluctuations vanish as E 3/2 , i.e., faster than the band energy near a node, rendering the semimetal phase stable.
The early field-theoretical point of view has been recently questioned by nonperturbative calculations [20,24], hinting that 3D gapless phases can become unstable due to the emergence of zero-energy states bound to statistically rare regions of the disorder potential landscape. According to this picture, the nodal DoS remains nonzero for arbitrarily weak disorder without any signature of singular behavior. Evidence for avoided quantum criticality (AQC) facilitated by localized nodal eigenstates has been provided by lattice simulations of a 3D Dirac model with uncorrelated on-site disorder [24]. Challenging these findings, Buchhold et al. noted that rare events are preceded by scattering resonances which always carry zero spectral weight at a node. Furthermore, as they are only possible for fine-tuned ("magical") impurity configurations, this would imply that the nodal DoS cannot be lifted at variance with the AQC scenario [27,28]. Their claim is backed by the exact solution of a Weyl node with a spherical impurity. These paradoxical findings have attracted significant attention recently [24,29,31,33] because they question the phase stability of incompressible 3D gapless phases in realistic conditions. Moreover, the driving mechanism for semimetal-to-metal transitions in the phase diagram of dirty 3D DWSMs remains elusive.
In this paper, we resolve this conundrum by tackling the spherical impurity problem using two complementary theoretical approaches. First, within a continuum model, we argue that an unforeseen nonanalytic behavior of scattering phase shifts at the node obstructs a direct use of Friedel's sum rule (FSR) [27,28]. This difficulty can be overcome by keeping track of the level statistics in systems of increasingly large volume at fixed impurity concentration. The puzzling behavior of the phase shifts is explained by the emergence and sudden disappearance of bound states upon tuning the impurity potential across a fine-tuned "magical value." Crucially, our level statistics analysis reveals that rare bound states are accompanied by a continuum of low-energy resonances surrounding the node in realistic material systems with a diversity of random short-range impurities. Such near-critical impurity configurations give effective statistical weight to magical impurities and ultimately endow the node with a finite average DoS. Second, we carry out high-precision numerical calculations in a lattice version of the problem, hosting one or more impurities. Remarkably, our real-space simulations not only unambiguously demonstrate the destabilization of a 3D DWSM by a diversity of random near-critical impurities but also quantitatively agree with the continuum theory predictions in the dilute impurity regime. These findings provide the missing link between continuum and lattice approaches to the DWSM theoretical puzzle and unambiguously pinpoint the newly unveiled statistical significance of near-critical impurities as the driving mechanism for AQC. Hence dilute impurities can only destabilize a DWSM provided their random parameters are drawn with a probability density which is nonzero on (at least) one magical value. Lastly, we note that subtleties in disordered Dirac systems have a long history [34][35][36][37][38]. For instance, in 2D d-wave superconductors, there are four lowenergy quasiparticle Dirac valleys, and scalar impurities are pair breaking. The latter induce resonances that, in the strong scattering limit, turn into sharp DoS peaks at E = 0 (Majorana zero modes) [39,40].
The remainder of the paper is organized as follows: In Sec. II, we present the theoretical tools for calculating the DoS correction induced by dilute spherical scalar impurities hosted within a single-node continuum model of a Dirac semimetal. We further argue our conclusions to remain valid in 3D Weyl semimetals. In Sec. III, we highlight the main caveats implied by a direct use of FSR and show how to obtain consistently the thermodynamic limit DoS change due to a finite (albeit small) concentration of impurities. In Sec. IV this theory is used to predict the conditions on which AQC holds in a Dirac semimetal. Our predictions are validated by high-resolution real-space calculations in Sec. V. Finally, in Sec. VI we summarize our main findings and highlight future directions for further study.

II. CONTINUUM THEORY
We start by reviewing the low-energy description of a noninteracting single-valley 3D DWSM. The Hamiltonian can be written as H 0 = vα · p, withh ≡ 1, α i = σ x ⊗ σ i , v being the Fermi velocity, and p =−ı∇ being the momentum operator. Here, σ i (i = x, y, z) denote Pauli matrices acting on internal spin space. Introducing a scalar impurity potential in the Hamiltonian breaks translation invariance, but if U (r) = U (|r|), rotational symmetry around the impurity center is preserved. For concreteness, we consider a spherical well or plateau potential, U (r) = λ (b −|r|) [20]. We note that this model is suitable to describe realistic multiple-valley Dirac or Weyl semimetals insofar as the impurity radius b is much larger than the lattice spacing (thus effectively suppressing intervalley scattering). The coupling of distinct Weyl sectors at each valley is also absent due to the scalar structure of the impurity potential. The eigenstates of H = H 0 + U (|r|) can be written as where j ∈{ 1 /2, 3 /2,...} and j z ∈{−j, − j + 1,..., j} are the total angular momentum quantum numbers, while κ =±1labels the eigenvalues of K = γ 0 · (2L · S − 1), i.e., κ ( j + 1 /2). Furthermore, ± j, j z (r) are orthonormal spin-1 /2 spherical harmonics, and f κ j (r)/g κ j (r) are radial functions. For nonzero energies, the latter are radial spherical waves with phase shifts introduced by the central potential (see Appendix A for additional details). In each j sector, the scattering phase shifts δ j induced by a spherical well or plateau are obtained by constraining the spinor to be continuous at r = b. One obtains, after a somewhat lengthy calculation, where Bessel functions of the first (second) kind. We underline that Eq. (2) is equivalent to that obtained in Refs. [20,28]f o r the Weyl equation. This is unsurprising because our 4 × 4 Dirac model is gapless and the impurity potential has a scalar structure (i.e., distinct Weyl sectors at each valley remain decoupled). However, Eq. (2) only defines δ j (ε, u)modπ . The ambiguity corresponds, at most, to a global change in the sign of the wave function. In order to obtain a unique definition of δ j (ε, u), one needs to choose a reference point: i.e., as the potential is switched off (u → 0), we require that the phase shifts vanish across the entire spectrum. A way to guarantee this is to enforce that δ j (ε →±∞, u) = −u [41][42][43], which is achieved by a trick explained in Appendix C. Since the analysis is qualitatively similar in all j sectors, in what follows we focus on the δ 1/2 (ε) phase shift (see Fig. 1) using the previous convention (for completeness, plots for other j's and around different magical u's are provided in Appendix C). For u = u j c , the phase shift is shown to have a physical π discontinuity at ε = 0, which marks the occurrence of zero-energy bound states [42]. For the 3D massless Dirac equation, bound states at ε = 0 can appear, for particular wells or plateaus, whenever a decoupling of the radial equations for f ± j (r > b)/g ± (r > b) occurs [20]. In this case, the admissible (asymptotically decreasing) solutions are simple power laws, As shown in Appendix A, such spinors are only continuous at r = b, if the potential satisfies J j (|u|) = 0. Hence zero-energy states are allowed in a single-impurity Dirac problem provided the parameters are fine-tuned, i.e., |λb/v|=u j c is a root of J j (x). The critical parameters {u j c } are dubbed magical values, as they would correspond to rare regions in a disordered landscape where nonperturbative zero-energy modes are possible [20,24]. Note that these are true (squared-normalizable) impurity bound states within the Dirac continuum and they have a degeneracy of 2(2 j + 1). These bound states manifest themselves as a π discontinuity at ε = 0 in the phase shifts when the parameter u crosses a critical value of that angular momentum channel (inset of Fig. 1). This is in accord with Levinson's theorem for Dirac particles [41,42], which states that given an appropriate convention, the number of bound states is encapsulated in discontinuous π jumps of the phase shifts at zero momentum.

III. IMPURITY-INDUCED CHANGE IN THE DoS
The change in the DoS induced by an isolated impurity is conventionally calculated using FSR, which measures the variation in the number of states per unit energy (an extensive quantity, not to be confused with ν per unit volume, hereafter denoted by ρ). Strikingly, the phase-shift discontinuity caused by the impurity bound state in the 3D DWSM problem precludes the direct use of FSR, a peculiar effect that has gone unnoticed in earlier studies. Therefore, to determine ν, we adopt a strategy based on counting states within a finite energy window adapted from Friedel's original reasoning [44]. First, we restrict the Dirac fermions to lie inside a finite sphere of radius R, S R .T h e Hermiticity of H is guaranteed if the Hilbert space is restricted to states with vanishing current across the spherical surface, ∂S R , that is, to a subspace where any two spinors and satisfyÕ ∂S R dS · [ † μ (r)α μν ν (r)] = 0. In Eq. (1), this is true if cos δ j (ε, u)J j (|ε|R) − sin δ j (ε, u)Y j (|ε|R) = 0. For each angular momentum sector, this condition quantizes the allowed energy levels, which we denote by ε j n .T h e number of ( j-sector) levels inside the energy window [ε 0 − ε/2,ε 0 + ε/2] is changed by the impurity due to an inwards or outwards migration of levels from regions of width ≃ δ j (ε ± ε /2, u)/R [up to O(R −2 )] near the respective boundaries. This mechanism is illustrated in Fig. 2(a). The variation in the number of j states inside the probing window is For a finite ε,E q .( 4) is accurate in the asymptotic limit R ≫ 1 and for ε 0 ± ε /2 = 0 as explained and illustrated in Appendix B. Next, we consider the intensive DoS ( ρ) induced by a finite density (c) of impurities in a volume V focusing on single scattering events. The neglect of quantumcoherent scattering by multiple impurities is justified in the dilute regime, where quantum interference corrections are suppressed by a factor of 1/(k F l ) [ 45], where l ∝ c −1 is the mean free path and k F = E /(vh) is the Fermi wavevector. This is an important assumption confirmed precisely by our numerical simulations below. Formally, the DoS is obtained by the limiting procedure, where i indexes the impurity and N (ε 0 , ε,{u i }, V )i st h ev a r i a t i o n in the total number of states. Assuming that {u i } are drawn from a probability density function p(u), the thermodynamic limit then reads where N j (ε 0 , ε,u) is given by Eq. (4). The final expression for the DoS variation due to a dilute set of random impurities, The order of limits in Eqs. (5) and (6)i se s s e n t i a l .T h e integration over u must be done prior to taking the ε → 0 + limit. This is reminiscent of lattice simulations, where the resolution parameter must be sent to zero only after the thermodynamic limit has been taken (see below). If δ j (ε, u) is differentiable at ε = ε 0 ,t h e ε → 0 + limit can be safely brought inside the integral, and one obtains ρ j (ε 0 ) = c(4 j + 2) ∂δ j (ε, u)/∂ε| ε=ε 0 u /π , where f u = du p(u) f (u), i.e., the familiar FSR. A direct application of FSR [Eq. (3)] was employed in Refs. [27,28] to determine the DoS in Weyl systems with statistical fluctuations of u around a critical value u c , leading to the conclusion that ρ j (0) = 0. This was inferred from the fact that ∂δ j (ε, u)/∂ε| ε=0 = 0 The curves were obtained using Eq. (3) with the phase shifts represented in Fig. 1. As one approaches u c , ρ 1/2 (ε) takes the form of a low-lying peak, which gets narrower and closer to ε = 0 (roughly conserving the area between zeros). The inset shows that ρ 1/2 (ε, u) is always equal to zero at ε = 0 for any noncritical u. (c) Depiction of N 1/2 (ε = 0,π)/ ε converging towards a distribution 4δ(u − π ) as ε → 0 + . (d) Theoretical prediction for ρ 1/2 (ε) due to dilute spherical impurities with a Gaussian diversity of width σ around u 1/2 c = π . The inset depicts the total DoS for 10 −6 b −1 impurities per volume.
for any u = u j c ; see Fig. 2(b). Since critical configurations (u = u j c ) have zero statistical measure, FSR would seemingly imply a vanishing average DoS at ε = 0. In the remainder of this paper, we show that the difficulty arising from the discontinuous δ j (ε, u = u c ) can be overcome by carefully accounting for the level statistics in the infinite volume limit of a DWSM with random impurities. Such a procedure does not alter the fate of ρ(ε = 0) in the presence of a finite concentration impurities fixed u [20,28]. However, the conclusions on the stability of a DWSM are changed dramatically if a continuous statistical diversity of "near-critical" impurities exists.

IV. NEAR-CRITICAL IMPURITIES LIFT THE NODAL DENSITY OF STATES
A finite concentration of exactly critical wells or plateaus would introduce a macroscopic number of nodal bound states. However, for a diversity of random impurities with potential strengths drawn from a probability distribution function p(u), such fine-tuned configurations appear with zero probability and cannot yield statistically significant contributions to the bulk nodal DoS [27,28]. Nevertheless, we find that the abundance of low-energy resonances due to near-critical configurations (u ≈ u j c ) provide such a contribution. The phase shifts of such impurities signal the emergence of the zero-energy bound states by a sharp resonance, namely, a quick π variation of δ j (ε)a su → u j c originated in the valence band, which moves towards ε = 0 and becomes sharper while always keeping δ j (0, u) = 0. At u = u c , the situation is delicate because δ j (ε) is no longer differentiable at ε = 0. In that case, one must work with Eq. (6) directly, and since there is a zero-energy π discontinuity in the phase shifts, a Dirac-δ distribution around the u j c,n emerges as the limit of plateau functions with a conserved integral equal to 4 (the degeneracy of the j = 1/2 single-impurity bound states). This limit is depicted in Fig. 2(a). An immediate implication is that a DWSM is unstable to dilute random impurities provided p(u j c ) = 0 for at least one critical u j c . Such a condition implies that the stability of a 3D semimetallic phase ultimately depends on the type of impurity model and the disorder statistics, i.e., whether it supports the resonant mechanism driven by a continuous distribution of near-critical impurities. These findings, supported below by accurate lattice simulations, show that dirty 3D DWSMs with near-critical impurities are inherently unstable, which sheds light on the previously reported AQC [24][25][26][29][30][31].
In Fig. 2(d), we plot the change in the j = 1/2D o S due to a dilute diversity of "near-critical impurities." The diversity is characterized by a Gaussian distribution p(u) = exp [−(u − π )/2σ 2 ]/ √ 2πσ around u c = π .T h eD o Si s clearly lifted around ε = 0, forming a sharp symmetrical bump. For this diversity model, the peak is Gaussian shaped near its center, and the corresponding area is conserved as σ → 0. In this limit, a 4δ(ε) distribution forms, i.e., all impurities are critical, each having a fourfold degenerate zeroenergy bound state.

V. LATTICE SIMULATIONS
Our prediction for the lifting of the DoS due to near-critical impurities has been based on a continuum model for a singlenode Dirac semimetal. However, real Dirac materials and numerical simulations live in the realm of lattice models, featuring several nodes and warped band structures. To validate our previous conclusions, we perform real-space simulations on a simple cubic lattice (L C of parameter a and linear size L) with a four-orbital Hamiltonian derived from the continuum Hamiltonian H [23][24][25], namely,  [47]. The energy resolution reads as η = π E /2M, where E is the bandwidth of the Hamiltonian matrix and M is the truncation order of the polynomial expansion [48,49]. The DoS is where ··· denotes disorder averaging and D = 4L 3 is the Hilbert space dimension. In order to simulate systems with a vanishing mean-level spacing, thereby performing calculations bounded only by η, we randomly sample over twisted boundaries [30,31]. This approach allowed the DoS to be calculated with unprecedented working spectral resolutions as low as η ≃ 4 × 10 −4 v/a, whose full convergence requires M ≈ 30 000 polynomials. More technical details are provided in Appendix D. Figure 3(a) shows the average DoS induced by critical impurities, ρ(E ) = ρ imp (E , u = π ) − 2E 2 /π 2 , in the dilute regime. The numerical data are compared with our analytical results [Eq. (6) and discussion thereafter], including the eightfold valley degeneracy, and properly convoluted with Gaussian functions of width η to mimic the finite numerical spectral resolution. The lifting of the DoS at the node and the underlying near-critical impurity mechanism are borne out by the spectral calculations, which show excellent quantitative agreement with the continuum theory, provided the impurity radius is large enough [see Fig. 3(a) and additional numerical evidence in Appendix D]. In Fig. 3(b), we present an analogous calculation for a system having several impurities inside the supercell. The impurities are placed randomly without superpositions, and their strengths drawn from a Gaussian distribution N (μ = π v/a, σ = 0.3 v/b). The continuum theory prediction for the low-energy bump in the DoS is reproduced in the diluted limit [see inset of Fig. 3(b)], and the law ρ(E = 0) ∝ c remains accurate up to 10 −6 impurities per unit cell. The overshooting for higher concentrations is due to multi-impurity effects, which become more effective as impurities are pushed closer together.

VI. CONCLUSIONS AND OUTLOOK
Here, we have shown that AQC must occur in 3D Dirac semimetals having dilute short-range scalar impurities, if their random parameters have a nonzero probability density at so-called magical values, where nodal bound states appear. These results were based on a continuum formulation of the problem treated at the single-impurity level and quantitatively confirmed by high-resolution lattice simulations in a gapless multivalley Dirac model hosting ≈10 −9 -10 −6 random scalar impurities impurities per unit cell. The perfect agreement between theory and numerical simulations gives confidence that the newly unveiled resonant mechanism stemming from diverse near-critical impurities is a crucial piece in the DWSM quantum criticality puzzle. Moreover, disparities with previous work [27,28] are explained by the presence of physical π jumps in the scattering phase shifts that prevent a direct use of FSR for diverse impurities around the aforementioned magical parameters. Similar conclusions are expected to hold for Weyl semimetals, as scalar impurities do not couple the different Weyl sectors in the infinite volume limit of our Dirac model. Meanwhile, our lattice nodal DoS calculations show a crossover from a dilute regime at very low impurity concentration (with DoS scaling linearly with c) to an intermediate impurity density regime (c 10 −6 impurities per unit cell), where the DoS diverges from the analytical prediction. This behavior can be traced to quantum-coherent multiple-impurity scattering events, which are neglected in our continuum theory.
A related, but nontrivial, question concerns the validity of these conclusions when dealing with lattice models having uncorrelated on-site disorder. In light of our theory, as well as earlier work [24,29], one reasonably expects the semimetallic phase to be unstable for unbounded distributions. However, such systems with highly concentrated and atomic-sized (onsite) impurities are exactly in the regime where the lattice results deviate from continuum predictions, hinting at yet further subtleties when relating the fate of disordered DWSM phases with rare-event bound states. Here, we provide technical details on the calculations leading from the Dirac eigenvalue problem in the presence of a single central scalar potential: H = H 0 + U (|r|). These results form the theoretical foundation for the main results presented in this paper. Since the methods employed are scattered around the existing literature [20,41,42], we provide a detailed description of procedures to make the presentation self-contained.

Derivation of the radial Dirac equations and radial eigenstates
The eigenvalue problem for an independent Dirac particle in the presence of a central potential corresponds to finding the solutions of where the repeated Greek indices are summed over the fourspinor components of the single-particle Dirac wave function. In particular, we are interested in the special case of a potential well or plateau, such that U (|r| < b) = λ and U (|r| b) = λ.
The first technical step towards solving Eq. (A1)i st o use a spherical coordinate system, (r,θ,ϕ), and achieve a separation of variables. The way to do this is well known in the relativistic quantum mechanics literature and is based on identifying the orbital and spin angular momentum operators for this system, which read where the matrices act in the Dirac spinor indices. These quantities are not conserved by the Hamiltonian H; however, we can build three mutually commuting observables out of L and S, which are conserved and uniquely define the spinor and angular structure of the eigenfunctions of H. These are and also K = γ 0 · (2L i S i − 1), which explicitly reads It is easy to verify that all three operators in Eqs. (A3a)-(A3c) commute among themselves and with H. Crucial for the latter is the fact that U (r) = U (r), which guarantees that the impurity does not break rotational symmetry around its center. Therefore a common eigenbasis of |J| 2 , J z , and K can be built and labeled by the set of quantum numbers j ∈{ 1 /2, 3 /2,...}, j z ∈{−j, − j + 1,..., j}, and κ =±1. The quantum number κ appears by solving K's eigenvalue problem, Using the previous operators, a general form for the eigenspinors indexed by the set ( j, j z ,κ)is where f κ j (r)/g κ j (r) are radial functions and κ (θ,ϕ) are spin- which in this form are orthonormalized in the unit sphere, i.e., π 0 sin θ dθˆ2 Besides the orthonormality condition of Eq. (A7), κ j, j z ( ) have some further useful properties, namely, σ ·r κ (θ,ϕ) = (σ ·r) 2 κ (θ,ϕ) = −κ (θ,ϕ), (A8a) where σ = (σ x ,σ y ,σ z ) is a vector of Pauli matrices and the scalar products are to be understood as a summation over spacial indices. Finally, we can proceed and write the Hamiltonian H explicitly as a differential operator in spherical coordinates. That way, it reads Using this form for H, one can plug spinors as in Eq. (A5)into the eigenvalue problem of Eq. (A1) and arrive at the following coupled systems of ordinary differential equations: In the case of the spherical well or plateau that concerns us, Eq. (A10) reduces to either outside of it. In Eqs. (A11) and (A12), we use dimensionless scales, namely, x = r/b, ε = Eb/v, and u = λb/v.Thesolutions inside the impurity (as long as ε = u) always have the general form where A ± are complex adjustable constants. Outside the impurity and for nonzero energy, one has instead where the choice of parametrization in the linear combination was made for convenience. Note that the exterior solutions feature both J n and Y n components, being always regular and physically admissible in their support (x 1 ) .N o w ,a l lw e must do is constrain the functions δ ± j (ε) such that the spinor (r) is continuous at the impurity's surface (x = 1). Using Eqs. (A14a)-(A14d), this implies that This equation is independent of κ, which allows us to define a unique function, δ j (ε, u), for both the κ =±sectors, which appears a single twofold degeneracy of the states in the problem. The previous facts justify Eq. (2) of Sec. II.
The previous analysis is valid for the entire spectrum, except at the important ε = 0 point. Here, the interior solutions are the same, but the radial system outside the impurity decouples, i.e., The latter allows for power-law solutions, of which the physically admissible ones (i.e., the ones decaying with x)a r eo f the form x j+1/2 . (A17b) Both these solutions, being joined continuously to the interior solutions of Eqs. (A13a)-(A13d) require that J j (|u|) = 0. This condition gives rise to a discrete set of parameters u,for which these bound-state solutions are allowed.
Finally, we remark that all the eigenstates determined here (the unbound and bound ones) have an intrinsic 2 j + 1 degeneracy due to the rotational invariance of the Hamiltonian. This degeneracy factor is explicitly taken into account in all calculations done in the main text.

APPENDIX B: SELF-ADJOINT RESTRICTION OF THE DIRAC HAMILTONIAN TO A FINITE SPHERE
To derive the relation between the scattering phase shifts and the change in the DoS due to a dilute diversity of impurities, we make explicit use of the restriction of H to a finite sphere of radius R ≫ b, i.e., S R . Restricting a continuum Hamiltonian to a finite volume of space generally makes its action on the original Hilbert space non-Hermitian. The way around this is to impose appropriate boundary conditions which restrict the original basis to a subset, generating a subspace inside of which the Hamiltonian preserves its Hermiticity. This is called taking a self-adjoint extension of H to a finite domain.
In the case of the Dirac Hamiltonian with a scalar potential, H =−ιvα · ∇ + U (r), the Hermiticity condition is imposed by guaranteeing that for any two Dirac spinor states 1 (r) and 2 (r), the following condition holds: After some straightforward manipulation, this condition can be cast into the equivalent form wheren = (n x , n y , n z ) is an outwards unit vector normal to the spherical surface ∂S R . This is precisely the condition 013183-7 FIG. 4. Plot of ε 1/2 = ε 1/2 n+1 − ε 1/2 n , the nearest-level spacings for a spherical impurity with u = 3.1867 calculated from the numerically found solutions of the boundary condition [Eq. (B3)]. The data points are represented as R 2 × ( ε 1/2 − π R −1 ), such that a collapse of different values of R is achieved. This collapse indicates that presented in Sec. III. Unsurprisingly, Eq. (B2) is easily interpreted as guaranteeing that no net particle current crosses the boundary of S R , which expresses particle conservation implied by Hermiticity.
Meanwhile, since all α matrices are composed of offdiagonal 2 × 2 blocks, one can easily see that Eq. (B2)i s satisfied whenever we impose either the first or the last two components of the Dirac spinors to be zero at ∂S R . Considering spinors of the form given in Eq. (A5), such a condition translates into either f + The other two combinations cannot be satisfied, as the zeros of Bessel functions of different j's never coincide. For the purposes of this work, we chose the first of these conditions (although the specific self-adjoint extension should not be relevant for any thermodynamic limit results). Finally, by using the general form of the exterior scattering solutions found earlier [Eqs. (A14a)-(A14d)], we arrive at our final form for the boundary condition, where R is measured in units of b.

Level spacing of central impurity Dirac Hamiltonian
First, we remark on an important consequence of the boundary condition in Eq. (B3). This condition imposes a quantization of energy levels, turning the continuous spectrum into a discrete one with a density of levels that scales with R. Provided that we are looking at finite energies (ε = 0) and with |ε|R ≫ 1, Eq. (B3) can be taken in its asymptotic form, namely, In the absence of an impurity, we have δ j (ε, 0) = 0, and the mesh of energy levels allowed by the boundary conditions (inagiven j sector) is simply ε j n ≈ nπ R + sgn(n) π 2R ( j + 1/2), with n ∈ Z. This yields a mean-level spacing which is uniform across the spectrum and equal to π R −1 . In the presence of the impurity (which induces energy-dependent phase shifts), Eq. (B4) does not seem to have a simple solution. However, if R is large enough such that δ j (ε, u) is a slowly varying function across an energy interval of width π R −1 , then one can say that the allowed energy levels are roughly which gives a correction to the mean-level spacing relative to the case u = 0, which is simply with an analogous expression for n < 0. Hence we conclude that the correction to the mean-level spacing due to a single impurity is always ∝O(R −2 ), which is subleading relative to the original π/R spacing. This result is exemplified by a numerical solution of Eq. (B3)i nF i g .4 and justifies our arguments on the number of states migrating in or out of an energy interval given in Sec. III.

APPENDIX C: A CONSISTENT DEFINITION OF THE SCATTERING PHASE SHIFTS
In this Appendix, we use the spherical Dirac eigenstates found earlier to define the scattering phase shifts in a way that allows a direct relation to the impurity-induced change in the density of states. We recover early results which explain the crucial zero-energy π discontinuity observed in our calculations as due to the appearance of critical bound states in the transition between the valence and conduction band. This connects our results to Levinson's theorem applied to noninteracting Dirac particles.
The exterior scattering wave functions of the gapless Dirac equation with a spherical well or plateau (of strength λ and radius b) were found to be of the form where N ± are complex normalization constants and δ j (E ,λb) are the energy-dependent scattering phase shifts. From the forms of Eqs. (C1) and (C2), it is clear that adding n × π (with integer n) to the phase shifts yields exactly the same spinor states, apart from irrelevant global sign changes. Meanwhile, the scattering phase shifts must always obey Eq. (A15), which guarantees continuity of the wave functions at r = b.H o wever, as explained in the main text, this condition does not uniquely define the functions δ j (ε, u), and a choice must be made concerning the reference situation relative to which the wave functions get dephased. A natural choice is to define δ j as the phase shift of the wave function relative to the case when u = 0. More precisely, we can think of a situation in which the central potential is adiabatically turned on and the instantaneous scattering eigenstates get progressively more dephased at all energies. This convention is known to be a useful one for Dirac fermions [41,42], and it can be achieved by enforcing δ j (ε →±∞, u) →−u. It is important to remark that this choice is needed for us to relate the change in the number of states inside a fixed spectral window with the phase shifts of scattering states in that window. In Fig. 5, we depict the lowestj phase shifts, δ 1/2 (ε, u) and δ 3/2 (ε, u), as a function of energy when u is close to a critical value. These curves were obtained using the previous convention for the phase shifts [δ j (ε →±∞, u) → −u], which can be guaranteed by the following numerical integration: Defining δ j by branches guarantees not only that the appropriate asymptotic convention is obeyed but also that the discontinuity due to zero-energy bound states is always avoided in the integrals. Finally, from Fig. 5 it is clear that a π discontinuity develops at ε = 0 when the impurity parameter is critical. This is the trademark of a zero-energy bound state since, according to Levinson's theorem for gapless Dirac particles, the number of bound states with well-defined j, j z , and κ is given as (see Ma and Ni [41]) which yields n j, j z ,κ=± (u = n j c ) = 0 and n j, j z ,κ=± (u = n j c ) = 1. This agrees with our earlier derivation of the zero-energy eigenstates in this system.

Technical description of the numerical method
Here, we provide some technical details on the numerical method used for calculating the density of states in the lattice model defined in Eq. (7) of Sec. V. As explained there, the calculations used a kernel polynomial method (KPM) [48], implemented in an efficient CPU parallelized framework developed by some of us (QUANTUM KITE [47]). We begin by outlining the basic elements of our numerical method.
Our aim is to calculate the intensive density of states (DoS) of a finite quantum lattice system with N degrees of freedom (in our case, N = 4L 3 , as we have a simple cubic lattice with side L and four orbitals per site). This quantity is given generically as ρ(ε)dε = 1 N α g α δ(ε − ε α )dε, where the summation is over eigenvalues of H and g α is the degeneracy of each level. For our numerical purposes, ρ(ε)dε is expanded in Chebyshev polynomials of the first kind, T n (x), yielding Finally, we remark that given a function f (x) which is approximated with a finite resolution in x,t h eK P Mapproximated function is described as the following convolution integral: This result is used for most of the analysis done on our realspace numerical results.

Lattice model and boundary conditions
In this brief section, we provide details and illustrate the lattice model used in all our numerical simulations. As referred to in the main text, our basic lattice Hamiltonian H D was obtained by discretizing the continuum Dirac Hamiltonian [with a scalar potential U (R)], H, in a simple cubic lattice with four orbitals per site. This tight-binding model Hamiltonian reads where a is the lattice parameter and † R = (c † R,A,↑ , c † R,A,↓ , c † R,B,↑ , c † R,B,↓ ) is a vector with on-site fermionic creation operators. Here, A and B stand for two different sublattices, while ↓ and ↑ are the two spin states in each orbital. Note that this convention for naming the local single-particle states is consistent with the previously defined intrinsic angular momentum operator S.I nF i g .6(a), we depict this real-space model in terms of its hoppings.
In the clean limit, U (R) = 0, this lattice model can be diagonalized by going to k space. After doing that, we obtain the particle-hole symmetric dispersion relation where both the conduction and valence bands are twofold degenerate. At half filling, this clearly reproduces a 3D Dirac semimetal, with eight valleys placed at the time-reversal invariant momenta (TRIM) of the first Brillouin zone. These are shown in Fig. 6(b). Near a TRIM, K D , the dispersion relation takes the form which is exactly the same as we had in our original continuum Hamiltonian. Nevertheless, the discretization of H 0 introduces a replication of the original fourfold degenerate Dirac cone into eight disconnected ones. Before ending this section, it is useful to calculate the normalized DoS of the clean lattice model, as it is used explicitly in the analysis of our numerical results. Using our previous definition, the intensive DoS for this system (assuming a simulated lattice with L 3 sites) reads v a sin 2 k x a + sin 2 k y a + sin 2 k z a . (D6) Due to particle-hole symmetry (ε →−ε), it suffices to evaluate the DoS at positive energies. Numerically, we can choose a regular mesh in the first Brillouin zone (FBZ) of the cubic lattice (equivalent to choosing a finite real-space cell) and determine the normalized DoS. This is shown in Fig. 7, together with the corresponding low-energy quadratic approximation. The latter is simply ρ 0 (E ) = 2E 2 /π 2 , which is the expression used in the main text.

Additional results for the resonances of a single sphere
In this section, we present additional details on the numerical results presented in the main text, together with additional results supporting our conclusions. We begin by presenting the numerical results for the change in the density of states due to a single extended sphere in the center of a simulated supercell of size L 3 . Despite simulating a single sphere, we are actually sampling over random realizations of boundary phase twists: This well-known technique helps the convergence of the KPM calculations, by eliminating the mean-level spacing from the problem. This method considers the computational domain as a supercell that gets repeated in a periodic cubic superlattice. In the large-L limit, periodicity artifacts eventually die out, and fluctuations around boundary-averaged values scale as ∝L −3/2 . Figure 8 shows numerical results for the ρ imp (E , u)L 3 = [ρ imp (E , u) − E 2 /2π 2 ]L 3 for values of u around π , the first positive critical value for a bound state with j = 1/2. These results are compared with theoretical curves (dashed black lines), obtained from the result of FSR, convoluted with a Gaussian, to account for the finite spectral resolution (η)impliedbythe numerical method. Note that Eq. (D7) includes a factor of 8, which accounts for the eight Dirac valleys existing in our lattice model, as well as a 1/4 due to the four orbitals per site in our lattice model. The numerically calculated DoS is then normalized by the number of states, 4L 3 , and the clean system has ρ(E ) ≈ 2E 2 /π 2 for E ≈ 0. As can be seen from Fig. 8, the agreement with the curves obtained from the continuum theory is perfect for spheres of radius b > 16 a with a concentration smaller than 256 −3 a −3 down to energy resolutions of meV. For spheres of radius b = 8 a, one already observes deviations from the continuum theory curves in the form of energy shifts (see bottom panels in Fig. 8).
In Fig. 9, we represent analogous high-energy-resolution numerical results for u ≈ 4.493 ... corresponding to the first resonance associated with j = 3/2. In the plots, one can also observe the next resonance (with j = 5/2) approaching the Dirac node. One can see that a radius of 16 a is not sufficiently large to have a complete agreement between the numerical peaks and the continuum theoretical curves. In the lower panels, the calculation is repeated for a larger radius of the spherical impurity (b = 22 a), and a perfect agreement is then obtained for j = 3/2.
Finally, it is important to analyze directly the case when the single impurity inside the supercell is at a critical value. In this case, we argue that uncoupled zero-energy eigenstates FIG. 8. Top: Plots of the change in the density of states due to a single spherical impurity of strength u = 2.8, 3.0, 3.4 v/a and radius b = 16 a, inside a simulated supercell of volume 256 3 a 3 . The vertical widths of the numerical curves are 95% statistical error bars, with respect to the simultaneous sampling over random KPM vectors and boundary phase twists. The agreement with the resolution-corrected theoretical curves is perfect over the entire range of resolutions used. Bottom: Results for a single spherical impurity of strength u = 2.8, 3.0 v/a and radius b = 8 a, inside a simulated supercell of volume 256 3 a 3 . Finite-size effects due to the discretization of the spherical impurity in the lattice are now visible as a shift of the peak away from the node. exist for the configuration, contributing as 8δ(E )/L 3 to the DoS (contributions coming from different valleys, as well as the factor of 4 due to the normalization to the total number of orbitals, are included). One can never see such a contribution numerically using the previous procedure, but we can analyze its emergence as a function of the spectral resolution. More precisely, we must have which is compared with numerical results (for u = π )i n Fig. 3. The agreement is perfect.
To close this section, we remark that the main conclusions to be drawn from the previous single-impurity results are threefold: (1) The continuum theory describes the DoS peaks corresponding to resonances associated with dilute near-critical spherical impurities, provided that this peak is located near the Dirac node (where the continuum theory holds), the radius of the spheres is large enough, and the distance between spheres is sufficiently large. (2) Largerj resonances require the discretized spheres to be larger in order to reproduce the continuum theory results for the same energy resolutions. (3) Numerically, one can observe the emergence of a Dirac δ at zero energy when the dilute impurities are all at critical values.

Details on the simulation of the average DoS for a system of random impurities
Here, we provide details on the generation of the random distribution of nonoverlapping spheres in the lattice used to produce the numerical results of Fig. 3(b). In order to do this, we started by considering a simulated supercell (with twisted boundaries) with 512 3 unit cells (≈536 000 000 orbitals), which from the results of the previous single-impurity simulations is sufficient to reproduce accurately the single sphere ρ(E ) at low energies if spheres of radius 16 a are considered. Then, we generate the potential associated with a regular cubic lattice composed by the centers of such (discretized) spheres inside the simulated cell. This procedure is equivalent to subdividing the original supercell side by an integer number, generating a set of identical subcells. FIG. 10. Scheme of the procedure used to generate a configuration of multiple random spheres inside the simulated supercell.
Finally, each of the generated central points is randomly displaced in three-dimensional space, and a potential strength is randomly chosen for each impurity inside the supercell. This procedure guarantees that there are no superpositions in any sample, as one restricts the random displacement of the centers to keep it inside the corresponding subcell. A schematic is depicted in Fig. 10. Once the full potential landscape inside the supercell is created, the remaining numerical procedure is identical to what was previously described.