\section{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\,\citep{Armitage2018}.
Several types of gapless systems featuring Dirac or Weyl points in
three-dimensional (3D) momentum space have been realized\,\citep{Borisenko2014,Han-Jin2017,Xu613}.
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
({\small{}$\mathcal{T}$}) or inversion symmetries ({\small{}$\mathcal{P}$})\,\citep{Armitage2018}
--- or crystal symmetries in {\small{}$\mathcal{TP}$}-symmetric
Dirac systems\,\citep{Young2012,Yang2014}. Thus, any clean DWSM
is an incompressible electron gas with a quadratically vanishing density
of states (DoS). \textcolor{black}{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\,\citep{PhysRevB.90.121107,PhysRevB.92.045104}.
In addition, weakly interacting DWSMs can display strongly renormalized
Fermi velocities, but the nodes' integrity and topology are expected
to remain robust\,\citep{Interact_PhysRevB.94.241102,Interact_PhysRevB.97.161102,Interact_PhysRevB.98.241102,Interact_PhysRevLett.113.136402,Roy16,Roy17,Roy18,maciejko2014}. }

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\,\citep{Fradkin1986_2,Shindou2009,Goswami2011,Yang2014,Nandkishore2014,Syzranov2015,Sbierski2016,Pixley2016,Pixley2016_2,pixley2016_3,Syzranov2016,Buchhold2018,Buchhold2018_2,Ziegler2018,Goncalves2020,Wilson2020,Kobayashi2020}.
An early result by Fradkin predicted that Dirac nodes are stable in
{\small{}$d\!=\!2\!+\!\epsilon$} dimensions, below some critical
disorder strength\,\citep{Fradkin1986_2}. 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{\small{} $W$}. Since, near a node, the de Broglie wavelength,
{\small{}$\lambda\!=\!\hbar v/E$}, largely exceeds the disorder correlation
length, the central limit theorem applies and the fluctuations around
the average potential inside a volume {\small{}$\lambda^{d}$ }must
scale as {\small{}$\delta V\!\!\propto\!W\!\lambda^{\!{\scriptscriptstyle -d/2}}$}.
In{\small{} $d\!=\!3$}, the fluctuations vanish as{\small{} $E^{{\scriptscriptstyle 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\,\citep{Nandkishore2014, Pixley2016_2},
hinting that 3D gapless phases can become unstable due to the emergence
of\emph{ 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\,\citep{Pixley2016_2}.
Challenging these findings, Buchhold \textit{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 (\textit{``magical}'') impurity configurations, this
would imply that the nodal DoS cannot be lifted at variance with the
AQC scenario\,\citep{Buchhold2018,Buchhold2018_2}. Their claim is
backed by the exact solution of a Weyl node with a spherical impurity.\textcolor{blue}{{}
}\textcolor{black}{These paradoxical findings have attracted significant
attention recently\,\citep{Pixley2016_2,roy_universal_2016,Ziegler2018,Wilson2020}
because they question the phase stability of incompressible 3D gapless
phases in realistic conditions. Moreover,}\textit{\textcolor{black}{{}
}}\textcolor{black}{the driving mechanism for semimetal-to-metal transitions
in the phase diagram of dirty 3D DWSMs remains elusive.}

\textcolor{black}{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
non-analytic behavior of scattering phase shifts at the node obstructs
a direct use of Friedel's sum rule (FSR)} \,\textcolor{black}{\citep{Buchhold2018,Buchhold2018_2}.
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}\textit{\textcolor{black}{\emph{
upon tuning the impurity potential across a fine-tuned}}}\textit{\textcolor{black}{{}
``magical value}}\textcolor{black}{''. Crucially, our level statistics
analysis reveals that rare bound states are accompanied by a }\textcolor{black}{\emph{continuum
of low-energy resonances}} surrounding the nod\textcolor{black}{e
in realistic material systems with a }\textcolor{black}{\emph{diversity}}\textcolor{black}{{}
of random} short-range impurities.\textcolor{red}{{} }Such near-critical
impurity configurations give effective statistical weight to \textit{\textcolor{black}{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\textcolor{black}{{} 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 }\textit{\textcolor{black}{statistical significance of near-critical
impurities}}\textcolor{black}{{} 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}\textit{\textcolor{black}{{} magical}}\textcolor{black}{{}
value.} Lastly, we note that subtleties in disordered Dirac systems
have a long history\,\citep{Atkinson2000,Hirschfeld2002, Pepin2001, Adagideli2002, Altland2002}.\,For
instance, in 2D $d$-wave superconductors, there are four low-energy
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)\,\citep{Balatsky95, Adagideli99}.

\textcolor{black}{\vspace{-0.1cm}
 }

\textcolor{black}{The remainder of the paper is organized as follows:\,In
Sec.\,\ref{sec:Continuum-theory}, 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.\,\ref{sec:Impurity-induced-change-in},
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.\,\ref{sec:Near-critical-impurities-lift}
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.\,\ref{sec:Lattice-simulations}.
Finally, in Sec.\,\ref{sec:Outlook} we summarize our main findings
and highlight future directions for further study.}

\vspace{-0.3cm}


\section{\label{sec:Continuum-theory}Continuum theory}

\vspace{-0.2cm}

We start by reviewing the low-energy description o\textcolor{black}{f
a noninteracting single-valley 3D DWSM. The Hamiltonian can be written
as $\mathcal{H}_{0}\!\!=\!v\,\mathbf{\boldsymbol{\alpha}}\!\cdot\mathbf{p}$,
with $\hbar\!\equiv\!1$, $\alpha^{i}\!=\!\sigma^{x}\!\!\otimes\!\!\sigma^{i}$,
$v$ being the Fermi velocity, and $\mathbf{p}\!=\!-\imath\boldsymbol{\nabla}$
being the momentum operator. Here, $\sigma_{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 $\mathcal{U}\!\left(\mathbf{r}\right)\!\!=\!\mathcal{U}\!\left(\abs{\mathbf{r}}\right)$,
rotational symmetry around the impurity center is preserved. For concreteness,
we consider a spherical well or plateau potential, $\mathcal{U}\!\left(\mathbf{r}\right)\!\!=\!\!\lambda\Theta\left(b\!-\!\left|\mathbf{r}\right|\right)$\,\citep{Nandkishore2014}.
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 $\mathcal{H}\!=\!\mathcal{H}_{0}\!+\!\mathcal{U}\!\left(\abs{\mathbf{r}}\right)$
can be written as}

\begin{equation}
\Psi_{j,j_{z}}^{\kappa}\!\left(\mathbf{r}\right)\!=\!\left[\frac{f_{j}^{\kappa}\!\left(r\right)}{r}\Theta_{j,j_{z}}^{-\kappa}\!\left(\hat{\mathbf{r}}\right),\frac{\iota g_{j}^{\kappa}\!\left(r\right)}{r}\Theta_{j,j_{z}}^{\kappa}\!\!\left(\hat{\mathbf{r}}\right)\right]^{\text{T}},\label{eq:SpinorAngularStructure-1}
\end{equation}
where $j\in\{\nicefrac{1}{2},\nicefrac{3}{2},\cdots\}$ and $j_{z}\in\left\{ -j,-j+1,\cdots,j\right\} $
are the total angular momentum quantum numbers, while $\kappa=\pm1$
labels the eigenvalues of $\mathcal{K}\!=\!\gamma^{0}\!\cdot\!(2\mathbf{L}\cdot\mathbf{S}\!-\!1)$,
i.e., $\kappa(j\!+\!\nicefrac{1}{2})$. Furthermore, $\Theta_{j,j_{z}}^{\pm}\!\!\left(\hat{\mathbf{r}}\right)$
are orthonormal spin-$\nicefrac{1}{2}$ spherical harmonics and $f_{j}^{\kappa}\left(r\right)$/$g_{j}^{\kappa}\left(r\right)$
are radial functions. For nonzero energies, the latter are radial
spherical waves with phase shifts introduced by the central potential
(see Appendix\,\eqref{sec:AppendixA} for additional details). In
each $j$ sector, the scattering phase shifts $\delta_{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, 
\begin{widetext}
\vspace{-0.4cm}

\begin{equation}
\tan\delta_{j}\left(\varepsilon,u\right)=\frac{\text{sgn}\left(\varepsilon\!-\!u\right)J_{j+1}\left(\abs{\varepsilon}\right)J_{j}\left(\abs{\varepsilon-u}\right)-\text{sgn}\left(\varepsilon\right)J_{j}\left(\abs{\varepsilon}\right)J_{j+1}\left(\abs{\varepsilon-u}\right)}{\text{sgn}\left(\varepsilon\right)\text{sgn}\left(\varepsilon-u\right)Y_{j+1}\left(\abs{\varepsilon}\right)J_{j}\left(\abs{\varepsilon-u}\right)-Y_{j}\left(\abs{\varepsilon}\right)J_{j+1}\left(\abs{\varepsilon\!-\!u}\right)},\label{eq:B/A-1-1-1}
\end{equation}

\vspace{-0.35cm}
 
\end{widetext}

\noindent where $u\!\equiv\!\lambda b/v_{\text{F}}$, $\varepsilon\!\equiv\!Eb/v_{\text{F}}\!\neq\!(0,u)$,
and $J_{n}(x)\,[Y_{n}(x)]$ are Bessel functions of the first (second)
kind. We underline that Eq.\,(\ref{eq:B/A-1-1-1}) is equivalent
to that obtained in Refs.\,\citep{Nandkishore2014,Buchhold2018_2}
for the Weyl equation. This is unsurprising because our $4\times4$
Dirac model is gapless and the impurity potential has a scalar structure
(i.e. distinct Weyl sectors at each valley remain decoupled). However,
Eq.\,\eqref{eq:B/A-1-1-1} only defines $\delta_{j}\!\left(\varepsilon,u\right)$\,mod$\pi$.
The ambiguity corresponds, at most, to a global change in the sign
of the wavefunction. In order to obtain a unique definition of $\delta_{j}\left(\varepsilon,u\right)$,
one needs to choose a reference point, i.e., as the potential is switched
off ($u\!\to\!0$), the require that the phase shifts vanish across
the entire spectrum\textcolor{black}{. A way to guarantee this is
to enforce that $\delta_{j}\left(\varepsilon\!\to\!\pm\infty,u\right)\!=\!-u$\,\citep{Ma85,Ma_2006}\,}\footnote{\textcolor{black}{As $\varepsilon\!\rightarrow\!\pm\infty$, one has
$\tan\delta_{j}\left(\varepsilon,u\right)\!\rightarrow\!-\tan u$.
This differs from the limiting behavior introduced in Ref.~\citep{Ma85},
$\lim\!{}_{\varepsilon\to\pm\infty}\!\delta_{j}\left(\varepsilon,u\right)\!\!=\!\!\mp u$.
The difference comes from the $\text{sgn}(\varepsilon)$ factors introduced
when defining the radial scattering wavefunctions {[}see Eqs.\,\eqref{eq:g+FreeSolution-3-2}---\eqref{eq:g-FreeSolutionOutside}
of Appendix\,\ref{sec:AppendixA}{]}}}\textcolor{black}{, which is achieved by a trick explained in Appendix\,\ref{sec:AppendixC}.
Since the analysis is qualitatively similar in all $j$ sectors, in
what follows we focus on the $\delta_{1/2}\!\left(\varepsilon\right)$
phase shift (see Fig.\,\ref{fig:PhaseShift - -Pi to Pi}) using the
previous convention (for completeness, plots for other $j$'s and
around different magical $u$'s are provided in Appendix\,\ref{sec:AppendixC}).
For $u\!\!=\!\!u_{c}^{j}$, the phase shift is shown to have a physical
$\pi$ discontinuity at $\varepsilon\!\!=\!\!0$, which marks the
occurrence of zero-energy bound states}\,\citep{Ma_2006}\textcolor{black}{.
For} the 3D massless Dirac equation, bound states at $\varepsilon\!=\!0$
can appear, for particular wells or plateaus, whenever a decoupling
of the radial equations for $f_{j}^{\pm}\left(r\!>\!b\right)/g^{\pm}\left(r\!>\!b\right)$
occurs\,\citep{Nandkishore2014}. In this case, the admissible (asymptotically
decreasing) solutions are simple power laws, $g_{j,k}^{+}\left(r\!>\!b\right)/f_{j,k}^{-}\left(r\!>\!b\right)\!=\!\mathcal{B}^{\pm}r^{-j-1/2}$
and $f_{j,k}^{+}\left(r\!>\!b\right)/g_{j,k}^{-}\left(r\!>\!b\right)\!\!=\!\!0$.
As shown in Appendix\,\eqref{sec:AppendixA}, such spinors are only
continuous at $r\!=\!b$, if the potential satisfies $J_{j}\left(\abs u\right)\!=\!0$.
Hence zero-energy states are allowed in a single-impurity Dirac problem
provided the parameters are fine-tuned, i.e., $\abs{\lambda b/v}\!=\!u_{c}^{j}$
is a root of $J_{j}\left(x\right)$. The critical parameters $\{u_{c}^{j}\}$
are dubbed \emph{magical values}, as they would correspond to rare
regions in a disordered landscape where nonperturbative zero-energy
modes are possible\,\citep{Nandkishore2014,Pixley2016_2}. Note that
these are true (squared-normalizable) impurity bound states within
the Dirac continuum and they have a degeneracy of $2\!\left(\!2j\!+\!1\!\right)$.
These bound-states manifest themselves as a $\pi$ discontinuity at
$\varepsilon\!=\!0$ in the phase shifts when the parameter $u$ crosses
a critical value of that angular momentum channel (inset of Fig.\,\ref{fig:PhaseShift - -Pi to Pi}).
This is in accord with Levinson's theorem for Dirac particles\,\citep{Ma85,Ma_2006}
which states that given an appropriate convention, the number of bound
states is encapsulated in discontinuous $\pi$ jumps of the phase
shifts at zero momentum.


\begin{figure}[t]
\vspace{-0.2cm}

\begin{centering}
\hspace{-0.1cm}\includegraphics[scale=0.165]{Fig1.pdf} 
\par\end{centering}
\vspace{-0.35cm}

\caption{\label{fig:PhaseShift - -Pi to Pi}Plot of the energy-dependent phase
shift $\delta_{{\scriptscriptstyle 1/2}}\!\left(\varepsilon,u\right)$
in accordance with the prescription $\delta_{j}\!\left(\varepsilon\!\to\!\pm\infty,u\right)\!=\!-u$.
Several values of $u$ are plotted around $u_{c}=\pi$\textcolor{black}{.\,The
inset depicts a close-up of the curves near $\varepsilon\!=\!0$.}\vspace{-0.8cm}
 }
\end{figure}

\vspace{-0.7cm}


\section{\label{sec:Impurity-induced-change-in}Impurity-induced change in
the DoS\textup{ }}

\vspace{-0.2cm}

\noindent The change in the DoS induced by an isolated impurity is
conventionally calculated using FSR,

\vspace{-0.55cm}

\begin{equation}
\Delta\nu(\varepsilon,u)=\frac{2}{\pi}\sum_{j=1/2}^{\infty}\left(2j+1\right)\frac{\partial\delta_{j}(\varepsilon,u)}{\partial\varepsilon},\label{eq:FSR}
\end{equation}

\textcolor{black}{\vspace{-0.2cm}
 }

\noindent \textcolor{black}{which measures the variation in the number
of states per unit energy (an extensive quantity, not to be confused
with $\Delta\nu$ per unit volume hereafter denoted by $\Delta\rho$).
Strik}ingly, 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 $\Delta\nu$, we adopt a strategy based on counting states
within a finite energy window adapted from Friedel's original reasoning\,\citep{Friedel1952}.
First, we restrict the Dirac fermions to lie inside a finite sphere
of radius $R$, $S_{R}$. The Hermiticity of $\mathcal{H}$ is guaranteed
if the Hilbert space is restricted to states with vanishing current
across the spherical surface, $\partial S_{R}$, that is, to a subspace
where any two spinors $\Psi$ and $\text{\ensuremath{\Phi}}$, satisfy
$\varoiint_{\partial S_{R}}\!d\mathbf{S}\!\cdot\!\left[\!\Psi_{\mu}^{\dagger}\left(\mathbf{r}\right)\!\boldsymbol{\alpha}_{\mu\nu}\Phi_{\nu}\left(\mathbf{r}\right)\!\right]\!=\!0.$
In Eq.\,\eqref{eq:SpinorAngularStructure-1}, this is true if $\cos\!\delta_{j}\left(\varepsilon,u\right)J_{j}\left(\abs{\varepsilon}\!R\right)\!-\!\sin\!\delta_{j}\left(\varepsilon,u\right)Y_{j}\left(\abs{\varepsilon}\!R\right)\!=\!0.$
For each angular momentum sector, this condition quantizes the allowed
energy levels, which we denote by $\varepsilon_{n}^{j}$. The number
of ($j$-sector) levels inside the energy window $[\varepsilon_{0}-\Delta\varepsilon/2,\varepsilon_{0}+\Delta\varepsilon/2]$
is changed by the impurity due to an inwards or outwards migration
of levels from regions of width $\simeq\delta_{j}\left(\varepsilon\pm\nicefrac{\Delta\varepsilon}{2},u\right)/R$
{[}up to $\mathcal{O}\left(R^{-2}\right)${]} near the respective
boundaries. This mechanism is illustrated in Fig. \ref{fig:LowEnergyDoS_Altland}(a).
The variation in the number of $j$ states inside the probing window
is

\vspace{-0.45cm}

\begin{align}
\Delta N_{j}\!\left({\scriptstyle \varepsilon_{0},\Delta\varepsilon,u}\!\right) & \!=\!\frac{2\left(2j\!+\!1\right)}{\pi}\!\left[\delta_{j}\!\left({\scriptstyle \varepsilon_{0}+\frac{\Delta\varepsilon}{2},u}\right)\!-\!\delta_{j}\left({\scriptstyle \varepsilon_{0}-\frac{\Delta\varepsilon}{2},u}\right)\right]\!.\label{eq:Number of States Variation}
\end{align}

\vspace{-0.2cm}

\noindent For a finite $\Delta\varepsilon$, Eq.\,\eqref{eq:Number of States Variation}
is accurate in the asymptotic limit $R\gg1$ and for $\varepsilon_{0}\pm\nicefrac{\Delta\varepsilon}{2}\neq0$
as explained and illustrated in Appendix\,\eqref{sec:AppendixB}.

Next, we consider the intensi\textcolor{black}{ve DoS ($\Delta\rho$)
indu}ced by a finite\textcolor{black}{{} density} ($c$) of impurities
in a volume $V$ focusing on \textcolor{black}{single scattering events.
The neglect of quantum-coherent scattering by multiple impurities
is justified in the dilute regime, where quantum interference corrections
are suppressed by a factor of $1/(k_{F}l)$\,\citep{altland_condensed_2010},
where $l\!\propto\!c^{-1}$ is the mean free path and $k_{F}=E/v\hbar$
is the Fermi wavevector. This is an important assumption confirmed
precisely by our numerical simulations below.}\textcolor{red}{{} }Formally,
the DoS is obtained by the limiting procedure, $\Delta\rho\left(\varepsilon_{0}\right)\!=\!\lim_{\Delta\varepsilon\to0^{+}}\lim_{V\rightarrow\infty}\left[\Delta N\left(\varepsilon_{0},\Delta\varepsilon,\left\{ u_{i}\right\} ,V\right)/V\Delta\varepsilon\right],$
where $i$ indexes the impurity and $\Delta N\left(\varepsilon_{0},\Delta\varepsilon,\left\{ u_{i}\right\} ,V\right)$
is the variation in the total number of states. Assuming that $\left\{ u_{i}\right\} $
are drawn from a probability density function $p\left(u\right)$,
the thermodynamic limit then reads

\vspace{-0.45cm}

\begin{align}
\sum_{i}\frac{\Delta N_{j}(\varepsilon_{0},\Delta\varepsilon,u_{i})}{V}\!\!\underset{{\scriptscriptstyle V\to\infty}}{\longrightarrow}\!c\!\int\!du\,p(u)\Delta N_{j}(\varepsilon_{0},\Delta\varepsilon,u),\label{eq:ThermodynamicLimitDiversity}
\end{align}
\vspace{-0.45cm}

\noindent where $\Delta N_{j}\!\left(\varepsilon_{0},\Delta\varepsilon,u\right)$
is given by Eq.~(\ref{eq:Number of States Variation}). The final
expression for the DoS variation due to a dilute set of random impurities,
$\Delta\rho\left(\varepsilon_{0}\right)\!=\!\!\sum_{j}\Delta\rho_{j}\left(\varepsilon_{0}\right)$,
is obtained from

\vspace{-0.45cm}

\begin{multline}
\Delta\rho_{j}(\varepsilon_{0})\!=\!c\!\lim_{{\scriptscriptstyle \Delta\varepsilon\to0^{+}}}\!\!\int du\,p\left(u\right)\frac{\Delta N_{j}\left(\varepsilon_{0},\Delta\varepsilon,u\right)}{\Delta\varepsilon}.\label{eq:Number of States Variation-2}
\end{multline}

\vspace{-0.3cm}

\noindent

\noindent 
\begin{figure}[t]
\vspace{-0.3cm}

\begin{centering}
\hspace{-0.25cm}\includegraphics[scale=0.235]{Fig2.pdf}
\par\end{centering}
\vspace{-0.35cm}

\caption{\label{fig:LowEnergyDoS_Altland}(a) Motion of energy levels triggered
by the central impurity. (b) Plots of $\Delta\rho_{{\scriptscriptstyle 1/2}}\!\left(\varepsilon,u\right)$
for selected values of $u$ around $u_{\text{c}}\!=\!\pi$. The curves
were obtained using Eq.(\ref{eq:FSR}) with the phase shifts represented
in Fig.\ref{fig:PhaseShift - -Pi to Pi}. As one approaches $u_{\text{c}}$,
$\Delta\rho_{1/2}\!\left(\varepsilon\right)$ takes the form of a
low-lying peak, which gets narrower and closer to $\varepsilon\!=\!0$
(roughly conserving the area between zeros). The inset shows that
$\Delta\rho_{1/2}\!\left(\varepsilon,u\right)$ is always equal to
zero at $\varepsilon\!=\!0$ for any non-critical $u$. (c) Depiction
of $\Delta N_{1/2}\!\left(\varepsilon\!=\!0,\pi\right)/\Delta\varepsilon$
converging towards a distribution $4\delta\!\left(u\!-\!\pi\right)$
as $\Delta\varepsilon\!\to\!0^{+}$. (d) Theoretical prediction for
$\Delta\rho_{1/2}\!\left(\varepsilon\right)$ due to dilute spherical
impurities with a Gaussian diversity of width $\sigma$ around $u_{c}^{{\scriptscriptstyle \!1\!/\!2}}\!\!=\!\pi$.
The inset depicts the total DoS for $10^{-6}b^{-1}$ impurities per
volume. \vspace{-0.8cm}
}
\end{figure}

\noindent The order of limits in Eqs. (\ref{eq:ThermodynamicLimitDiversity})
and (\ref{eq:Number of States Variation-2}) is essential. The integration
over $u$ must be done prior to taking the $\Delta\varepsilon\!\to\!0^{{\scriptscriptstyle +}}$
limit\textcolor{black}{. 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 $\delta_{j}\!\left(\varepsilon,u\right)$
is differentiable at $\varepsilon\!=\!\varepsilon_{0}$, the $\Delta\varepsilon\!\to\!0^{{\scriptscriptstyle +}}$
limit can be safely brought inside the integral, and one obtains $\Delta\rho_{j}\left(\varepsilon_{0}\right)\!=\!c\left(4j\!+\!2\right)\!\av{\partial\delta_{j}\left(\varepsilon,u\right)\!/\!\partial\varepsilon|_{\varepsilon=\varepsilon_{0}}\!}_{\!u}\!/\pi,$
where $\langle f\rangle_{u}\!=\!\int\!du\,p\left(u\right)\!f(u)$,
i.e., the familiar FSR. A direct application of FSR\textcolor{black}{\,{[}Eq.\,\eqref{eq:FSR}{]}}
was employed in Refs.\,\citep{Buchhold2018,Buchhold2018_2} to determine
the DoS in Weyl systems with statistical fluctuations of $u$ around
a critical value $u_{c}$, leading to the conclusion that $\Delta\rho_{j}\!\!\left(0\right)\!=\!0$.
This was inferred from the fact that $\partial\delta_{j}\left(\varepsilon,u\right)/\partial\varepsilon|_{\varepsilon=0}\!=\!0$
for any $u\!\neq\!u_{c}^{j}$; see Fig.\,\ref{fig:LowEnergyDoS_Altland}(b).
\textcolor{black}{Since critical configurations ($u\!=\!u_{c}^{j})$
have zero statistical measure, FSR would seemingly imply a vanishing
average DoS at $\varepsilon\!=\!0$. In the remainder of this paper,
we show that the difficulty arising from the discontinuous $\delta_{j}\!\left(\!\varepsilon,\!u\!=\!u_{\text{c}}\!\right)$
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 $\rho\left(\varepsilon\!=\!0\right)$
in the presence of a finite concentration of impurities with a fixed
$u$\,\citep{Buchhold2018_2,Nandkishore2014}. However, the conclusions
on the stability of a DWSM are changed dramatically if a }\textit{\textcolor{black}{continuous
statistical diversity of ``near-critical'' impurities}}\textcolor{black}{{}
exists. }

\vspace{-0.6cm}


\section{\label{sec:Near-critical-impurities-lift}Near-critical impurities
lift the nodal density of states}

\vspace{-0.2cm}

\noindent \textcolor{black}{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\,\citep{Buchhold2018,Buchhold2018_2}.
Nevertheless, we find that low-energy resonances due to }\textit{\textcolor{black}{near-critical}}\textcolor{black}{{}
configurations ($u\!\approx\!u_{c}^{j}$) 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 $\pi$ variation
of $\delta_{j}\!\left(\varepsilon\right)$ as $u\!\to\!u_{c}^{j}$
originated in the valence band, which moves towards $\varepsilon\!=\!0$
and becomes sharper while always keeping $\delta_{j}\!\left(0,u\right)\!=\!0$.\,At
$u\!=\!u_{c}$, the situation is delicate because $\delta_{j}\!\left(\varepsilon\right)$
is no longer differentiable at $\varepsilon\!=\!0$.\,In that case,
one must work with Eq.\,\eqref{eq:Number of States Variation-2}
directly, and since there is a zero-energy $\pi$ discontinuity in
the phase shifts, a Dirac-$\delta$ distribution around the $u_{\text{c},n}^{j}$
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.\,\ref{fig:LowEnergyDoS_Altland}(a).\,An
immediate implication is that a DWSM is unstable to dilute random
impurities provided $p(u_{c}^{j})\!\neq\!0$ for at least one critical
$u_{c}^{j}$.\,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\,\citep{Pixley2016_2,pixley2016_3,Syzranov_2018,Ziegler2018,Goncalves2020,Wilson2020}.}

\textcolor{black}{In Fig.\,\ref{fig:LowEnergyDoS_Altland}(d), we
plot the change in the $j\!\!=\!\!1/2$ DoS due to a dilute diversity
of}\textit{\textcolor{black}{{} ``near-critical impurities}}\textcolor{black}{''.
The diversity is characterized by a Gaussian distribution $p(u)\!=\!\exp\left[-\left(u\!-\!\pi\right)/2\sigma^{2}\right]/\sqrt{2\pi}\sigma$
around $u_{\text{c}}\!=\!\pi$. The DoS is clearly lifted around $\varepsilon\!=\!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 $\sigma\!\to\!0$. In this limit, a $4\delta\!\left(\varepsilon\right)$
distribution forms, i.e., all impurities are critical, each having
a fourfold degenerate zero-energy bound state.}

\vspace{-0.5cm}


\section{\label{sec:Lattice-simulations}Lattice simulations}

\vspace{-0.2cm}

Our prediction for the lifting of the DoS due to near-critical impurities
has been based on a continuum model for a single-node 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 ($\mathcal{L}_{C}$ of parameter $a$ and
linear size $L$) with a four-orbital Hamiltonian derived from the
continuum Hamiltonian $\mathcal{H}$\,\citep{Pixley2016,Pixley2016_2,pixley2016_3},
namely,

\begin{figure}[t]
\begin{centering}
\hspace{-0.3cm}\includegraphics[scale=0.24]{Fig3.pdf}
\par\end{centering}
\vspace{-0.3cm}

\caption{\label{fig:LowEnergyDoS}(a) DoS change due to an impurity of critical
strength $u\!=\!\pi\,v/a$ and radius $16\,a$ inside a supercell
of $256^{3}$ sites. Vertical widths are $95\%$ statistical error
bars, and dashed lines are the continuum theory predictions. (b) Plot
of $\protect\av{\rho\!\left(E\!=\!0\right)}_{u}$ with several impurities
of radius $16\,a$ inside the simulated supercell of $512^{3}$ sites
for different resolutions $\eta$. The gray line is the dilute regime
prediction. The inset shows converged $\protect\av{\rho\!\left(E\right)}_{u}$
for three concentrations against the predictions of Fig.\,\ref{fig:LowEnergyDoS_Altland}(d)
(black lines).\vspace{-0.7cm}
}
\end{figure}


\vspace{-0.5cm}

\begin{align}
\negthickspace\!H & =\!\negthickspace\sum_{{\scriptscriptstyle \mathbf{R}\in\mathcal{L_{\text{C}}}}}\negthickspace\left[\frac{\imath v}{2a}\Psi_{\mathbf{{\scriptscriptstyle R}}}^{\dagger}\!\cdot\!\alpha^{j}\!\cdot\!\Psi_{{\scriptscriptstyle \mathbf{R}+a\hat{e}_{j}}}\!+\!\frac{\mathcal{U}\left(\mathbf{R}\right)}{2}\Psi_{\mathbf{{\scriptscriptstyle R}}}^{\dagger}\!\cdot\!\Psi_{\mathbf{{\scriptscriptstyle R}}}\!+\!\text{H.c.}\right]\!.\label{eq:LatticeModel}
\end{align}

\vspace{-0.2cm}

\noindent The DoS is calculated by means of accurate Chebyshev polynomial
expansions of the resolvent operator~$\delta(E\!-\!H)$\,\footnote{\textcolor{black}{In all numerical calculations, the variable $E$
is to be understood as a dimensionless energy, measured in units of
$a/\hbar v$.}} in very large systems as implemented in the open-source quantum transport
code $\mathtt{QUANTUM}\;\mathtt{KITE}$\,\citep{Joao2020}. The energy
resolution reads as $\eta\!=\!\pi\Delta E/2M$, where $\Delta E$
is the bandwidth of the Hamiltonian matrix and $M$ is the truncation
order of the polynomial expansion\,\citep{Weise2006,Ferreira2015}.
The DoS is obtained from $\rho(E)\!=\!\lim_{\eta\rightarrow0}\lim_{D\rightarrow\infty}D^{-1}\textrm{Tr}[\langle\delta_{\eta}(E\!-\!H)\rangle]$,
where $\langle...\rangle$ denotes disorder averaging and $D\!=\!4L^{3}$
is the Hilbert space dimension. \textcolor{black}{In order to simulate
systems with a vanishing mean-level spacing, thereby performing calculations
bounded only by $\eta$, we randomly sample over twisted boundaries\,\citep{Goncalves2020,Wilson2020}.
}This approach allowed the DoS to be calculated with unprecedented
working spectral resolutions as low as $\eta\!\simeq\!4\!\times\!10^{-4}v/a$\textcolor{black}{,
whose full convergence requires $M\!\approx\!30\,000$ polynomials.
}More technical details are provided in Appendix\,\ref{sec:AppendixD}.

\textcolor{black}{Figure\,\ref{fig:LowEnergyDoS}(a) shows the average
DoS induced by critical impurities, $\Delta\rho\left(E\right)\!=\!\rho_{\text{imp}}\left(E,u\!=\!\pi\right)\!-\!2E^{2}/\pi^{2}$,
in the dilute regime. The numerical data are compared with our analytical
results {[}Eq.\,\eqref{eq:Number of States Variation-2} and discussion
thereafter{]}, including the eightfold valley degeneracy, and properly
convoluted with Gaussian functions of width $\eta$ 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.\,\ref{fig:LowEnergyDoS}(a) and additional numerical
evidence in Appendix\,\ref{sec:AppendixD}{]}. In Fig.\,\ref{fig:LowEnergyDoS}(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 $\text{\ensuremath{\mathcal{N}}}\!\!\left(\mu\!=\!\pi\,v/a\!\!,\sigma\!=\!0.3v/b\right)$.
The continuum theory prediction for the low-energy bump in the DoS
is reproduced in the diluted limit {[}see inset of Fig.\,\ref{fig:LowEnergyDoS}(b){]},
and the law $\rho\!\left(E\!=\!0\right)\!\propto\!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. }

\vspace{-0.3cm}


\section{\label{sec:Outlook}Conclusions and Outlook}

\vspace{-0.2cm}

\textcolor{black}{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 }\textit{\textcolor{black}{magical}}\textcolor{black}{{}
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 }\textcolor{black}{\emph{multivalley}}\textcolor{black}{{}
Dirac model hosting $\approx10^{-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\,\citep{Buchhold2018,Buchhold2018_2}
are explained by the presence of physical $\pi$ jumps in the scattering
phase shifts that prevent a direct use of FSR for diverse impurities
around the aforementioned }\textit{\textcolor{black}{magical}}\textcolor{black}{{}
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\gtrsim10^{-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.}

\textcolor{black}{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\,\citep{Pixley2016_2,Ziegler2018}, one reasonably expects
the semimetallic phase to be unstable for unbounded distributions.
However, such systems with highly concentrated and atomic-sized (on-site)
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.}

\vspace{-0.3cm}


\section{\textup{Acknowledgments}}

\vspace{-0.2cm}

We acknowledge support by the Portuguese Foundation for Science and
Technology through Strategic Funding No. UIDB/04650/2020, Projects
No. POCI-01-0145-FEDER-028887 (J.P.S.P. and J.M.V.P.L) and No. CEECIND/02936/2017
(B.A.), and Grant No. PD/BD/142774/2018 (J.P.S.P.). A.F. acknowledges
financial support from the Royal Society through a Royal Society University
Research Fellowship. The numerical calculations were performed on
the Viking Cluster, a high performance computing\,(HPC) facility
provided by the University of York. We are grateful for support from
the University of York's HPC service and their team. The authors are
grateful to A. Altland, M. Buchhold, M. Gonçalves, S.\,M. João, A.
Antunes, and J.\,M.\,B. Lopes dos Santos for valuable discussions.
The authors finally thank the anonymous referees for useful comments
which helped to improve the clarity of this work.

\vspace{-0.25cm}


