# Nodal vacancy bound states and resonances in three-dimensional Weyl
semimetals

J. P. Santos Pires<EMAIL_ADDRESS>Departamento de Física e Astronomia,
Faculdade de Ciências da Universidade do Porto, Rua do Campo Alegre, s/n,
4169-007 Porto, Portugal  Centro de Física das Universidades do Minho e do
Porto (CF-UM-UP) and Laboratory of Physics for Materials and Emergent
Technologies LaPMET, University of Porto, 4169-007 Porto, Portugal  S. M. João
Departamento de Física e Astronomia, Faculdade de Ciências da Universidade do
Porto, Rua do Campo Alegre, s/n, 4169-007 Porto, Portugal  Centro de Física
das Universidades do Minho e do Porto (CF-UM-UP) and Laboratory of Physics for
Materials and Emergent Technologies LaPMET, University of Porto, 4169-007
Porto, Portugal  Aires Ferreira  Department of Physics and York Centre for
Quantum Technologies, University of York, YO10 5DD, York, United Kingdom  B.
Amorim  Centro de Física das Universidades do Minho e do Porto (CF-UM-UP) and
Laboratory of Physics for Materials and Emergent Technologies LaPMET,
Universidade do Minho, 4710-057 Braga, Portugal  J. M. Viana Parente Lopes
<EMAIL_ADDRESS>Departamento de Física e Astronomia, Faculdade de Ciências da
Universidade do Porto, Rua do Campo Alegre, s/n, 4169-007 Porto, Portugal
Centro de Física das Universidades do Minho e do Porto (CF-UM-UP) and
Laboratory of Physics for Materials and Emergent Technologies LaPMET,
University of Porto, 4169-007 Porto, Portugal

###### Abstract

The electronic structure of a cubic $\mathcal{T}$-symmetric Weyl semimetal is
analyzed in the presence of atomic-sized vacancy defects. Isolated vacancies
are shown to generate nodal bound states with $r^{{\scriptscriptstyle-2}}$
asymptotic tails, even when immersed in a weakly disordered environment. These
states show up as a significantly enhanced nodal density of states which, as
the concentration of defects is increased, reshapes into a nodal peak that is
broadened by inter-vacancy hybridization into a comb of satellite resonances
at finite energies. Our results establish point defects as a crucial source of
elastic scattering that leads to nontrivial modifications in the electronic
structure of Weyl semimetals.

##  I Introduction

With the advent of three-dimensional (3D) topological insulators [1, 2], the
search for topological semimetals emerging at the transition between gapped
phases of matter has flourished. Rather than being a fine-tuned situation, it
was envisaged by Murakami [3] that, without inversion symmetry, topological
phase transitions can proceed through an intermediate stage, in which a pair
of two-fold degenerate band-crossing points moves around the first Brillouin
zone until it finally merges together and gives rise to a new gapped phase.
Such a stable gapless state was dubbed a Weyl semimetal (WSM) [4, 5] because
low-energy excitations around these band-crossings are described by a
decoupled pair of ($3\\!+\\!1)$–dimensional Weyl equations of opposite
chirality [6]. Later on, such a topological gapless phase was also shown to be
possible in centrossymetric crystals, so long as time-reversal symmetry is
broken [7, 8, 9, 10] (magnetic WSM). Crucially, in all cases, the band-
crossings form pairs of point-like sources (or sinks) of Berry flux in
momentum space, analogous to the well-known “diabolical points” described by
Berry [11] in a generic two-level quantum system. Therefore, isolated Weyl
nodes are topologically protected degeneracies in the electronic band
structure that are robust to parametric changes of the Hamiltonian.

The topological character of WSMs yields important physical consequences, from
the existence of surface Fermi arcs [7, 12, 13, 14, 15] that connect Weyl
nodes in the surface-projected first Brillouin zone (fBZ), to the remarkable
condensed matter realization of QED’s chiral anomaly [16, 17]. The latter
drives distinctive unconventional transport effects, such as a negative
longitudinal magnetoresistance [18, 19], a giant in-plane Hall effect [20, 21]
and the chiral magnetic effect [22]. Nonetheless, perhaps the most remarkable
property of a WSM is its resilience to the effects of unavoidable
perturbations, such as disorder or crystal defects. From a theoretical
standpoint, the study of disorder effects in both spectral [23, 24, 25, 26,
27, 28] and transport properties [29, 30, 31, 32, 33] of WSMs have been the
subject of intense research. A big focus was placed on the effects of random
potentials that can yield non-Anderson quantum criticality at a finite
disorder strength [29, 34, 35]. As the system is driven through this critical
point, the semi-metallic character of the nodal single-electron states gets
destroyed long before they become exponentially localized at the Anderson
transition [24].

In contrast to conventional Anderson transitions [36, 37], the disorder-
averaged nodal density of states (nDoS) in disordered WSMs is deemed an
appropriate order parameter by field-theoretical calculations [38, 39, 30,
40], as well as the numerical observation of its sharp power-law growth above
some critical disorder strength [41, 24, 25]. However, recent studies of
nonperturbative instantonic effects have revealed that rare disorder
configurations lift the nDoS and round-out its critical behavior [42, 43, 28],
thus challenging the conventional scenario. A physical picture was then put
forward by Nandkishore et al [42], who associated the nDoS lift to smooth
rare-regions of a random potential landscape that can sporadically bound
eigenstates at the nodal energy, giving way to an exponentially small but
_nonzero_ nDoS. Despite being a controversial proposal [25, 44, 45, 46, 26,
27], the _avoided_ quantum criticality due to rare events was eventually
confirmed in subsequent numerical studies [25, 26, 27] and it is now believed
to be the most general scenario [28].

In spite of these numerous theoretical studies regarding the effects of random
potential disorder, very little is known about the role played by point
defects and other common disorder sources. Currently, time-reversal symmetric
Weyl fermions can be realized as low-energy quasiparticles in a myriad of
materials, most notably within the TaAs cubic family (also including NbAs, TaP
and NbP) [47], which can be grown as single crystals using chemical vapor
transport techniques [48]. In the growth process, lattice defects are likely
to form [49, 50] and, as demonstrated in previous experimental studies based
on transmission electron microscopy [51] and Raman scattering [52], even high-
quality samples generally host a considerable density of defects, mostly
vacancies and stacking faults. Adding to their natural occurrence, vacancy
defects can also be artificially induced by means of particle irradiation [53,
54], a well-tested technique previously used to generate defects in graphene
[55] and two-dimensional semiconductors [56, 57]. Since point defects can
significantly change the electronic structure of materials, a study of their
impact as a source of disorder in WSMs opens up interesting possibilities.
Promising results were reported by Xing et al. [58], where atomic vacancies
hosted by the (magnetic) WSM $\text{Co}_{3}\text{Sn}_{2}\text{S}_{2}$ were
linked to the presence of exotic localized spin–orbit polaron states on its
surface. In this paper, we push this line forward by theoretically analyzing
the electronic properties of Weyl fermions in the presence of point defects.
More specifically, we characterize the electronic wavefunctions and
corresponding density of states (DoS) of a lattice $\mathcal{T}$-symmetric
Weyl semimetal with finite concentrations of randomly distributed atomic-sized
vacancies with one (half-vacancy) or two (full-vacancy) orbitals missing from
the defect sites.

The remainder of this paper is organized as follows. In Sec. II, we introduce
our working model and a projected Green’s function formalism (pGF) that is
used to calculate the vacancy-induced DoS deformation and show that
algebraically decaying nodal bound states appear for isolated half- and full-
vacancies. In Sec. III, the existence of nodal bound states is further
verified by Lanczos diagonalization [59, 60, 61] (LD) of lattices containing
an isolated vacancy. The robustness of these states to an additional weakly
disordered environment is also discussed. In Sec. IV, we analyze the averaged
DoS of a WSM with a finite concentration of vacancies, employing a combination
of LD and spectral methods [62, 63, 64]. While confirming that localized
eigenstates still appear and enhance the value of the DoS around the Weyl
node, our results further show that inter-vacancy hybridization quickly
broadens the nodal peak in the DoS, forming a comb of symmetrically-placed
subsidiary sharp resonances for a moderate concentration of defects. Finally,
Sec. V summarizes our key results and gives an outlook.

##  II Modeling an Isolated Vacancy in a Weyl Semimetal

A lattice vacancy is a common crystalline defect [65]. When a crystal is
formed some sites are not properly occupied by the corresponding atoms,
creating a proportion of vacant sites [66, 67] that act as a source of
disorder. In the language of tight-binding Hamiltonians, a vacancy can be
modeled by removing one or more Wannier orbitals from a randomly chosen
lattice site. We start by determining the effects of introducing a single
lattice half- or full-vacancy in a two-band model of a WSM. We employ a
particle-hole symmetric model that lives in a simple cubic lattice
($\mathcal{L}$) and features a low-energy dispersion relation with eight
isotropic Weyl nodes pinned to the time-reversal invariant momenta of the
cubic fBZ (see Fig. 1). The lattice Hamiltonian [28] may be written as

Figure 1: Band structure of the clean lattice WSM model along the
$\boldsymbol{k}$-space path indicated in the inset. The locations of the eight
non-equivalent Weyl cones are represented as well.

$\mathcal{H}_{0}\\!=\\!\frac{\hbar
v}{2a}\sum_{\mathbf{R}\in\mathcal{L}}\sum_{{\scriptscriptstyle
i=x\\!,\\!y\\!,\\!z}}\\!\\!\left[\Psi_{\mathbf{R}}^{\dagger}\\!\cdot\\!\sigma^{i}\\!\\!\cdot\\!\Psi_{\mathbf{R}+a\hat{\mathbf{x}}_{i}}\\!-\\!\text{h.c.}\right]$
(1)

where $a$ is the lattice parameter, $v$ is the Fermi velocity,
$\hat{\mathbf{x}}_{i}\\!=\\!(\hat{x},\hat{y},\hat{z})$ are cartesian unit
vectors, $\boldsymbol{\sigma}$ is the vector of $2\\!\times\\!2$ Pauli
matrices and
$\Psi_{\mathbf{R}}^{\dagger}\\!=\\![c_{\mathbf{R},1}^{\dagger},c_{\mathbf{R},2}^{\dagger}]$
is a local two-orbital fermionic creation operator. Equipped with this lattice
description, the vacancy defects are implemented in two distinct ways. In our
pGF calculations below, lattice vacancies are created by canceling all
hoppings at the defect site, which leaves behind uncoupled zero energy Wannier
states. In contrast, when the system is analyzed using spectral methods or LD
[Secs. III and IV], the Hilbert space’s dimension is effectively reduced by
iterating with vectors orthogonal to the removed orbitals.

###  II.1 Clean Lattice Green Function and Nodal Point Symmetries

Before diving into the analysis of the electronic structure of WSMs with
vacancy defects, we first establish some basic results. The Bloch Hamiltonian
of the clean system admits the simple representation

$\mathcal{H}_{0}(\mathbf{k})\\!=\\!i\hbar
v\boldsymbol{\sigma}\\!\cdot\\!\boldsymbol{\sin}\,a\mathbf{k},$ (2)

with $\boldsymbol{\sin}\,a\mathbf{k}\equiv\left(\sin ak_{x},\sin ak_{y},\sin
ak_{z}\right)$, and which yields the dispersion relation represented in Fig.
1. The clean retarded lattice Green’s function (LGF), defined formally as
$G_{0}(E,\mathbf{R}_{j}\\!-\\!\mathbf{R}_{i})\\!=\\!\left[E\\!+\\!i0^{{\scriptscriptstyle+}}\\!\\!\\!-\\!\mathcal{H}_{0}\right]_{\mathbf{R}_{j},\mathbf{R}_{i}}^{{\scriptscriptstyle-1}}$,
can be written, in terms of dimensionless quantities, as

$\boldsymbol{G}_{0}(\varepsilon,\boldsymbol{\Delta
R})\\!=\\!\\!\\!\int_{{\scriptscriptstyle\\![-\\!\pi,\pi]^{3}}}\\!\\!\\!\frac{d^{{\scriptscriptstyle\text{(\\!3\\!)}}}\boldsymbol{k}}{8\pi^{3}}\frac{\varepsilon\\!+\\!i\boldsymbol{\sigma}\\!\cdot\\!\boldsymbol{\sin
k}}{\varepsilon^{2}\\!+\\!\left|\boldsymbol{\sin
k}\right|^{2}}e^{-i\mathbf{k}\cdot\boldsymbol{\Delta R}},$ (3)

where $\varepsilon\\!=\\!Ea/\hbar v\\!+\\!i\eta$ is the dimensionless energy
(shifted by an imaginary amount $\eta$), $\boldsymbol{k}$ is the crystal
momentum (in units of $a^{{\scriptscriptstyle-1}}$), $\boldsymbol{\Delta
R}\\!=\\!\left(n_{x},n_{y},n_{z}\right)$, and $n_{i}\\!\in\\!\mathbb{Z}$ are
indices that label particular sites. Equation (3) can be expressed in terms of
four basic integrals over the domain $[-\pi,\pi]^{3}$, i.e.,

$\boldsymbol{G}_{0}(\varepsilon,\\!\boldsymbol{\Delta
R})\\!=\\!\mathcal{I}_{\varepsilon}^{0}\\!\\!\left(\boldsymbol{\Delta
R}\right)-\\!\\!\\!\\!\\!\sum_{j=x,y,z}\\!\\!\\!\\!i\sigma^{j}\,\mathcal{I}_{\varepsilon}^{j}\\!\left(\boldsymbol{\Delta
R}\right),$ (4)

where the complex-valued integrals are

Figure 2: Scheme of the local perturbation defined in Eq. (9).

$\displaystyle\mathcal{I}_{\varepsilon}^{0}\\!\left(\\!\boldsymbol{\Delta
R}\right)\\!$
$\displaystyle=\\!\\!\int_{{\scriptscriptstyle\\![-\\!\pi,\pi]^{3}}}\\!\\!\\!\\!d^{{\scriptscriptstyle\text{(\\!3\\!)}}}\\!\boldsymbol{k}\frac{\varepsilon\,e^{ik_{x}n_{x}}e^{ik_{y}n_{y}}e^{ik_{z}n_{z}}}{8\pi^{3}\left(\varepsilon^{2}+\left|\boldsymbol{\sin
k}\right|^{2}\right)}\,,$ (5a)
$\displaystyle\mathcal{I}_{\varepsilon}^{x}\\!\left(\\!\boldsymbol{\Delta
R}\right)\\!$
$\displaystyle=\\!\\!\\!\int_{{\scriptscriptstyle\\![-\\!\pi,\pi]^{3}}}\\!\\!\\!\\!\\!d^{{\scriptscriptstyle\text{(\\!3\\!)}}}\\!\boldsymbol{k}\frac{\sin\\!k_{x}e^{ik_{x}n_{x}}e^{ik_{y}n_{y}}e^{ik_{z}n_{z}}}{8\pi^{3}\left(\varepsilon^{2}+\left|\boldsymbol{\sin
k}\right|^{2}\right)},$ (5b)

with $\mathcal{I}_{\varepsilon}^{y,z}\\!\left(\\!\boldsymbol{\Delta R}\right)$
being obtained from
$\mathcal{I}_{\varepsilon}^{x}\\!\left(\\!\boldsymbol{\Delta R}\right)$ by a
cyclic permutation of the set $(n_{x},n_{y},n_{z})$. These integrals can be
numerically evaluated with high precision (see Appendix A) and display useful
symmetry properties. First, the parity of the integrand dictates that
$\boldsymbol{G}_{0}(E,\boldsymbol{\Delta R})\\!=\\!0$ whenever the set
$(n_{x},n_{y},n_{z})$ contains more than one odd integer. This small-scale
property of the LGF can be traced back to the existence of eight non-
equivalent valleys which are specific to this lattice model. Additionally,
there are non-spatial symmetries which can be deduced from Eqs. (5a)-(5b),
most notably

$\displaystyle\mathcal{I}_{\varepsilon}^{0}\\!\left(\\!\boldsymbol{\Delta
R}\right)\\!$
$\displaystyle=\\!\left[\mathcal{I}_{\varepsilon^{*}}^{0}\\!\left(\\!\boldsymbol{\Delta
R}\right)\right]^{*}$ (6a)
$\displaystyle\mathcal{I}_{\varepsilon}^{j}\\!\left(\\!\boldsymbol{\Delta
R}\right)\\!$
$\displaystyle=\\!-\\!\left[\mathcal{I}_{\varepsilon^{*}}^{j}\\!\left(\\!\boldsymbol{\Delta
R}\right)\right]^{*}$ (6b)

and also

$\text{Re}\\!\left[\mathcal{I}_{i0^{+}}^{0}\\!\left(\\!\boldsymbol{\Delta
R}\right)\right]\\!=\\!\\!\\!\\!-\lim_{\eta\to
0^{+}}\\!\int_{\text{[-{\hbox{\pi}},\leavevmode{p}]}^{3}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!d^{{\scriptscriptstyle\text{(\\!3\\!)}}}\boldsymbol{k}\frac{\eta\,\sin\left(\boldsymbol{k}\cdot\boldsymbol{\Delta
R}\right)}{8\pi^{3}\left(\eta^{2}\\!+\\!\left|\boldsymbol{\sin
k}\right|^{2}\right)}\\!=\\!0,$ (7)

by employing the cubic symmetry of the fBZ. Together, these three properties
imply that the LGF at $E\\!=\\!0$ can be represented in the simple form

$\boldsymbol{G}_{0}(0,\\!\boldsymbol{\Delta
R})\\!=\\!\\!\\!\\!\\!\\!\sum_{j=x,y,z}\\!\\!\\!\\!\sigma^{j}\text{Im}\\!\left[\mathcal{I}_{i0^{+}}^{j}\\!\\!\left(\boldsymbol{\Delta
R}\right)\right],$ (8)

which is non-zero if and only if $\boldsymbol{\Delta
R}\\!=\\!(n_{x},n_{y},n_{z})$ features a single odd integer.

###  II.2 Projected Green’s Function for a Lattice Vacancy

Within a lattice description, a vacancy can be modeled by removing hoppings
connecting one (or several) orbitals within a unit cell to its neighbors. With
no loss of generality, let us consider a vacancy at the origin,
$\boldsymbol{R}\\!=\\!\boldsymbol{0}$. The Hamiltonian is then
$\mathcal{H}\\!=\\!\mathcal{H}_{0}\\!+\\!\mathcal{V}$ with

$\displaystyle\mathcal{V}\\!$ $\displaystyle=\\!-\frac{\hbar
v}{2a}\\!\\!\sum_{{\scriptscriptstyle
i=x\\!,\\!y\\!,\\!z}}\\!\\!\left[\Psi_{\boldsymbol{0}}^{\dagger}\\!\cdot\\!\sigma^{i}\\!\\!\cdot\\!\Psi_{\boldsymbol{0}+a\hat{\mathbf{x}}_{i}}\\!\\!-\\!\Psi_{\boldsymbol{0}}^{\dagger}\\!\cdot\\!\sigma^{i}\\!\\!\cdot\\!\Psi_{\boldsymbol{0}-a\hat{\mathbf{x}}_{i}}\\!\\!-\\!\text{h.c.}\right].$
(9)

This operator has the advantage of having a finite support, i.e., it acts only
on sites $\Omega\\!=\\!\left\\{\mathbf{0},\pm a\hat{\mathbf{x}}_{1},\pm
a\hat{\mathbf{x}}_{2},\pm a\hat{\mathbf{x}}_{3}\right\\}$ that form the
octahedron shown in Fig. 2. Such local perturbations to a lattice model can be
conveniently studied by using the pGF method. Treating $\mathcal{V}$ as a
perturbation, we obtain a set of Dyson’s equations for the system’s Green’s
function, $\boldsymbol{G}(\varepsilon)$, in the presence of the vacancy, i.e.,

$\displaystyle\boldsymbol{G}(\varepsilon)$
$\displaystyle\\!=\\!\boldsymbol{G}_{0}(\varepsilon)\\!+\\!\boldsymbol{G}(\varepsilon)\\!\cdot\\!\mathcal{V}\\!\cdot\\!\boldsymbol{G}_{0}(\varepsilon)$
(10a) $\displaystyle\boldsymbol{G}(\varepsilon)$
$\displaystyle\\!=\\!\boldsymbol{G}_{0}(\varepsilon)\\!+\\!\boldsymbol{G}_{0}(\varepsilon)\\!\cdot\\!\mathcal{V}\\!\cdot\\!\boldsymbol{G}(\varepsilon),$
(10b)

where $\cdot$ denotes the matrix product defined in the full Hilbert space. To
solve these equations, we proceed in two steps: i) By projecting them into
$\Omega$, we can solve for those entries of the full
$\boldsymbol{G}(\varepsilon)$, i.e.,

$\overline{\boldsymbol{G}}(\varepsilon)\\!=\\!\left[\mathcal{I}\\!-\\!\overline{\boldsymbol{G}_{0}}(\varepsilon)\\!\cdot\\!\mathcal{V}\right]^{{\scriptscriptstyle-1}}\\!\\!\\!\cdot\overline{\boldsymbol{G}_{0}}(\varepsilon),$
(11)

where $\overline{\boldsymbol{G}_{0}}(\varepsilon)$ is the clean lattice
Green’s function restricted to $\Omega$, with
$\mathcal{I}\\!-\\!\overline{\boldsymbol{G}_{0}}(\varepsilon)\\!\cdot\\!\mathcal{V}$
defined within the (finite-dimensional) Hilbert subspace of $\Omega$, and ii)
The continuation of $\boldsymbol{G}(\varepsilon)$ to the exterior of $\Omega$
is obtained via

$\boldsymbol{G}(\varepsilon)\\!=\\!\boldsymbol{G}_{0}(\varepsilon)\\!+\\!\boldsymbol{G}_{0}(\varepsilon)\\!\cdot\\!\overline{\mathcal{T}_{\varepsilon}}\\!\cdot\\!\boldsymbol{G}_{0}(\varepsilon),$
(12)

where

$\displaystyle\overline{\mathcal{T}_{\varepsilon}}\\!$
$\displaystyle=\\!\mathcal{V}\\!+\\!\mathcal{V}\\!\cdot\\!\overline{\boldsymbol{G}}(\varepsilon)\\!\cdot\\!\mathcal{V}$
(13)
$\displaystyle=\mathcal{V}\\!\cdot\\!\left[\mathcal{I}\\!-\\!\overline{\boldsymbol{G}}_{0}(\varepsilon)\\!\cdot\\!\mathcal{V}\right]^{{\scriptscriptstyle-1}}$

is the projected $T$-matrix of the vacancy.

The pGF method provides access to the electronic structure of an isolated
impurity or defect embedded in an otherwise perfect infinite crystal [68, 45].
In the present section, we are interested in i) the DoS change induced by a
vacancy, and ii) the possibility that a WSM can host nodal bound states around
the vacancy. The emergence of zero-energy modes due to disorder is not obvious
given the absence of non-spatial symmetries in our model; note that
$\mathcal{H}_{0}$ belongs to the orthogonal Wigner-Dyson class (class AI in
the Altland-Zirnbauer ten-fold classification [28]). This is to be contrasted
to the well-studied case of graphene (chiral orthogonal BDI class [69, 70]),
which supports zero-energy states localized around point defects whose
peculiar spectral and transport properties have been linked to the underlying
chiral symmetry of that model [63, 71].

One of the most readily available observable from the pGF formalism is the
change in the (extensive) DoS due to the vacancy, $\delta\nu(\varepsilon)$.
This can be evaluated by means of equation,

$\displaystyle\delta\nu(\varepsilon)\\!$ $\displaystyle=\\!\frac{1}{\pi
n_{b}}\text{Im}\\!\left(\text{Tr}\left[\boldsymbol{G}_{0}(\varepsilon)\\!-\\!\boldsymbol{G}(\varepsilon)\right]\right)$
$\displaystyle=\\!-\frac{1}{\pi
n_{b}}\text{Im}\\!\left(\text{Tr}\left[\mathcal{V}\\!\cdot\\!\left[\mathcal{I}\\!-\\!\boldsymbol{G}_{0}(\varepsilon)\\!\cdot\\!\mathcal{V}\right]^{{\scriptscriptstyle-1}}\\!\\!\\!\cdot\left(\boldsymbol{G}_{0}(\varepsilon)\right)^{2}\right]\right)$
$\displaystyle=\\!\frac{1}{\pi
n_{b}}\text{Im}\\!\left(\text{tr}\left[\overline{\mathcal{T}_{\varepsilon}}\cdot\frac{d}{d\varepsilon}\boldsymbol{G}_{0}(\varepsilon)\right]\right),$
(14)

where $\overline{\mathcal{T}_{\varepsilon}}$ is the projected $T$-matrix,
$\text{Tr}\left[\cdots\right]$ is the trace operation over all degrees of
freedom, $\text{tr}\left[\cdots\right]$ is a trace over the support of
$\mathcal{V}$, $\Omega$, and $n_{b}$ is the number of orbitals per unit cell
($n_{b}\\!\\!=\\!\\!2$ in the WSM model). In the last step to obtain Eq. (14),
we have used Eq. (13) as well as the identity
$d\boldsymbol{G}_{0}(\varepsilon)/d\varepsilon\\!=\\!-\left(\boldsymbol{G}_{0}(\varepsilon)\right)^{2}$.

Next, we discuss briefly how to extract information on bound states within the
pGF framework. We start by writing the Lippmann-Schwinger equation for a
scattering state $\left|\Psi_{\varepsilon}\right\rangle$,

$\left|\Psi_{\varepsilon}\right\rangle\\!=\\!\left|\Psi_{\varepsilon}^{0}\right\rangle\\!+\\!\boldsymbol{G}_{0}(\varepsilon)\\!\cdot\\!\mathcal{V}\left|\Psi_{\varepsilon}\right\rangle,$
(15)

where $\left|\Psi_{\varepsilon}^{0}\right\rangle$ is an eigenstate of the
unperturbed system, which is the parent extended state of
$\left|\Psi_{\varepsilon}\right\rangle$. In contrast, a bound state can exist
without any parent eigenstate of the clean Hamiltonian. Thus, an eigenstate
bound by $\mathcal{V}$ at an energy $\varepsilon_{b}$ must be a solution of

$\left|\Psi_{\varepsilon_{b}}\right\rangle\\!=\\!\boldsymbol{G}_{0}(\varepsilon_{b})\\!\cdot\\!\mathcal{V}\left|\Psi_{\varepsilon_{b}}\right\rangle.$
(16)

Since the perturbation $\mathcal{V}$ has a finite support, one can once again
consider the projected version of Eq. (15),

$\left|\xi_{\varepsilon}\right\rangle\\!=\\!\overline{\boldsymbol{G}_{0}}(\varepsilon)\\!\cdot\\!\mathcal{V}\left|\xi_{\varepsilon}\right\rangle,$
(17)

where $\left|\xi_{\varepsilon_{b}}\right\rangle$ is the restriction of
$\left|\Psi_{\varepsilon_{b}}\right\rangle$ to the support $\Omega$. Thereby,
any bound state must obey the condition,

$\left[\mathcal{I}\\!-\\!\overline{\boldsymbol{G}_{0}}(\varepsilon_{b})\\!\cdot\\!\mathcal{V}\right]\left|\xi_{\varepsilon_{b}}^{b}\right\rangle\\!=\\!0,$
(18)

which means that its projected wavefunction must belong to the kernel of the
operator,
$\mathcal{I}\\!-\\!\overline{\boldsymbol{G}_{0}}(\varepsilon_{b})\\!\cdot\\!\mathcal{V}$.
Outside the support of $\mathcal{V}$ the wavefunction may be reconstructed
using

$\Psi_{\alpha}^{b}\\!\left(\mathbf{R}\right)\\!=\\!\left\langle\mathbf{R},\alpha\mid\Psi_{\varepsilon_{b}}^{b}\right\rangle\\!=\\!\left\langle\mathbf{R},\alpha\right|\\!\boldsymbol{G}_{0}(\varepsilon_{b})\\!\cdot\\!\mathcal{V}\left|\xi_{\varepsilon_{b}}^{b}\right\rangle,$
(19)

which may or may not amount to a normalizable state, depending on the
asymptotic behavior of the clean LGF. In this context, since we are looking
for zero-energy modes ($\varepsilon_{b}\\!=\\!0$), any state obeying Eq. (18)
is guaranteed to be square-normalizable in 3D space with an algebraic tail
$\propto r^{{\scriptscriptstyle-2}}$. The latter is the long-distance behavior
of $\boldsymbol{G}_{0}(0,\boldsymbol{\Delta R})$, as obtained in the continuum
limit.

####  II.2.1 Full-Vacancies in a Weyl Semimetal

We now apply the general formalism described above to the case of a full-
vacancy where both orbitals are removed from a particular lattice cell. To
perform the calculation, it comes in handy to order the sites of $\Omega$ as
$\left\\{\boldsymbol{0},a\hat{x},a\hat{y},a\hat{z},-a\hat{x},-a\hat{y},-a\hat{z}\right\\}$.
With this ordering, one obtains $\mathcal{V}$ as the matrix

$\mathcal{V}\\!=\\!\frac{\hbar
v}{2a}\left[\begin{array}[]{ccc}\mathbb{O}_{{\scriptscriptstyle
2\\!\times\\!2}}&\\!\\!\\!-\boldsymbol{\sigma}&\boldsymbol{\sigma}\\\
\boldsymbol{\sigma}^{T}&\mathbb{O}_{{\scriptscriptstyle
6\\!\times\\!6}}&\mathbb{O}_{{\scriptscriptstyle 6\\!\times\\!6}}\\\
\\!\\!\\!-\boldsymbol{\sigma}^{T}&\mathbb{O}_{{\scriptscriptstyle
6\\!\times\\!6}}&\mathbb{O}_{{\scriptscriptstyle
6\\!\times\\!6}}\end{array}\right],$ (20)

within the projected subspace, while, by exploiting the symmetries of
$\mathbf{G}_{0}$ [Eqs. 6a and 6b], we are able to further express the pGF as

$\overline{\boldsymbol{G}_{0}}(0)\\!=\\!\frac{a}{i\hbar
v}\\!\\!\left[\,\,\begin{array}[]{ccc}\mathbb{O}_{{\scriptscriptstyle
2\\!\times\\!2}}&g_{0}\boldsymbol{\sigma}&\\!\\!\\!\\!-g_{0}\boldsymbol{\sigma}\\\
\\!\\!\\!\\!-g_{0}\boldsymbol{\sigma}^{T}&\mathbb{O}_{{\scriptscriptstyle
6\\!\times\\!6}}&\mathbb{O}_{{\scriptscriptstyle 6\\!\times\\!6}}\\\
g_{0}\boldsymbol{\sigma}^{T}&\mathbb{O}_{{\scriptscriptstyle
6\\!\times\\!6}}&\mathbb{O}_{{\scriptscriptstyle
6\\!\times\\!6}}\end{array}\,\right],$ (21)

where $g_{0}\\!=\\!\mathcal{I}_{i0^{+}}^{x}\\!\left(1,0,0\right)$. These
simple matrices can then be used to build the operator
$\mathcal{I}\\!-\\!\overline{\boldsymbol{G}_{0}}(0)\\!\cdot\\!\mathcal{V}$,
whose determinant takes the remarkably simple form,

$\det\left(\mathcal{I}\\!-\\!\overline{\boldsymbol{G}_{0}}(0)\\!\cdot\\!\mathcal{V}\right)\\!=\\!\left(i\\!-\\!3g_{0}\right)^{4}.$
(22)

Equation (22) has a clear physical interpretation: a four-fold degenerate root
appears for $g_{0}\\!=\\!\nicefrac{{i}}{{3}}$, corresponding to an extra pair
of non-trivial bound states extending into the lattice. For consistency,
$\det\left(\mathcal{I}\\!-\\!\overline{\boldsymbol{G}_{0}}(0)\\!\cdot\\!\mathcal{V}\right)$
would have to display a two-fold degenerate root, corresponding to the
subspace of both orbitals to be removed from the lattice. Surprisingly, the
degeneracy appears doubled here, which indicates that two nodal bound states
must exist as a mathematical property of a full-vacancy in our lattice model.
In addition, a full diagonalization of
$\mathcal{I}\\!-\\!\overline{\boldsymbol{G}_{0}}(0)\\!\cdot\\!\mathcal{V}$
yields the following projected wavefunctions for these states:

$\displaystyle\left|\xi_{1}^{b}\right\rangle$
$\displaystyle=\frac{1}{\sqrt{6}}\left[\left|\hat{\mathbf{x}},1\right\rangle-i\left|\hat{\mathbf{y}},1\right\rangle-\left|\hat{\mathbf{z}},2\right\rangle\right.$
(23a)
$\displaystyle\qquad-\left.\left|-\hat{\mathbf{x}},1\right\rangle+i\left|-\hat{\mathbf{y}},1\right\rangle+\left|-\hat{\mathbf{z}},2\right\rangle\right],$

$\displaystyle\left|\xi_{2}^{b}\right\rangle$
$\displaystyle=\frac{1}{\sqrt{6}}\left[\left|\hat{\mathbf{x}},2\right\rangle\right.+i\left|\hat{\mathbf{y}},2\right\rangle+\left|\hat{\mathbf{z}},1\right\rangle$
(23b)
$\displaystyle\qquad-\left.\left|-\hat{\mathbf{x}},2\right\rangle-i\left|-\hat{\mathbf{y}},2\right\rangle-\left|-\hat{\mathbf{z}},1\right\rangle\right],$

where $\left|\mathbf{R},\alpha\right\rangle$ are the local Wannier states
(here, $\alpha$ indexes the orbital). Upon a reconstruction, these states have
the following real-space wavefunction outside $\Omega$:

$\displaystyle\Psi_{1}^{b}\\!\left(\mathbf{R}\right)\\!$
$\displaystyle=\\!\frac{it\sqrt{6}}{2}\left[\begin{array}[]{c}G_{0}^{12}\\!\left(0,\mathbf{R}\right)\\\
G_{0}^{22}\\!\left(0,\mathbf{R}\right)\end{array}\right]$ (24c)
$\displaystyle\Psi_{2}^{b}\\!\left(\mathbf{R}\right)\\!$
$\displaystyle=\\!\frac{it\sqrt{6}}{2}\left[\begin{array}[]{c}G_{0}^{11}\\!\left(0,\mathbf{R}\right)\\\
G_{0}^{21}\\!\left(0,\mathbf{R}\right)\end{array}\right]\\!,$ (24f)

where $G_{0}^{\alpha\beta}$ are spinor components of $\mathbf{G}_{0}$.

In a similar manner, we can determine the DoS change caused by a vacancy
defect. To do this, we require the clean pGF at all energies, which has the
following matrix structure,

$\overline{\boldsymbol{G}_{0}}(\varepsilon)\\!=\\!\frac{a}{i\hbar
v}\\!\\!\left[\,\,\begin{array}[]{ccc}f_{\varepsilon}\mathbb{I}_{{\scriptscriptstyle
2\\!\times\\!2}}&g_{\varepsilon}\boldsymbol{\sigma}&\\!\\!\\!\\!-g_{\varepsilon}\boldsymbol{\sigma}\\\
\\!\\!\\!\\!-g_{\varepsilon}\boldsymbol{\sigma}^{T}&\,\,f_{\varepsilon}\mathbb{I}_{{\scriptscriptstyle
6\\!\times\\!6}}&h_{\varepsilon}\mathbb{I}_{{\scriptscriptstyle
6\\!\times\\!6}}\\\
g_{\varepsilon}\boldsymbol{\sigma}^{T}&\,\,h_{\varepsilon}\mathbb{I}_{{\scriptscriptstyle
6\\!\times\\!6}}&f_{\varepsilon}\mathbb{I}_{{\scriptscriptstyle
6\\!\times\\!6}}\end{array}\,\right]\\!\\!,$ (25)

as imposed by the aforementioned symmetries of our WSM model. In Eq. (25),
$f_{\varepsilon}\\!=\\!\mathcal{I}_{\varepsilon}^{0}\left(0,0,0\right)$,
$g_{\varepsilon}\\!=\\!\mathcal{I}_{\varepsilon}^{x}\left(1,0,0\right)$ and
$h(\varepsilon)\\!=\\!\mathcal{I}_{\varepsilon}^{0}\left(2,0,0\right)$ are
dimensionless functions of the energy variable. Having the pGF, we can now
build the $T$-matrix of an isolated full-vacancy and employ Eq. (14) to
obtain,

$\\!\\!\\!\delta\nu(\varepsilon)\\!=\\!\frac{3a}{\pi\hbar
v}\text{Im}\\!\\!\left[\\!\frac{f_{\varepsilon}\left(h_{\varepsilon}^{\prime}\\!-\\!2f_{\varepsilon}^{\prime}\right)\\!+\\!f_{\varepsilon}^{\prime}h_{\varepsilon}\\!-\\!4g_{\varepsilon}^{\prime}\left(i\\!+\\!3g_{\varepsilon}\right)}{3f_{\varepsilon}^{2}\\!-\\!3f_{\varepsilon}h_{\varepsilon}\\!+\\!2\left(i\\!+\\!3g_{\varepsilon}\right)^{2}}\\!\right]\\!\\!.$
(26)

The functions $f_{\varepsilon}$, $g_{\varepsilon}$ and $h_{\varepsilon}$ and
their derivatives were calculated numerically (see Appendix A) and the
resulting DoS is shown in Fig. 3. As expected, a single (full) vacancy causes
a negative correction to the DoS across the entire band, which is consistent
with an overall transfer of spectral weight to the emergent bound states
around the vacant site. The integral of this curve is exactly $-2$, as the
number of continuum states is not conserved, i.e., two states (per orbital)
appear as a vacancy bound-states and other two are removed from the Hilbert
space.

####  II.2.2 Half-Vacancies in Weyl Semimetals

We now consider the effects of an isolated half-vacancy in which only one
orbital is removed from each cell. If the orbital $1$ is vacant, the
perturbation reads

$\displaystyle\mathcal{V}_{\\!1}\\!$ $\displaystyle=\\!-\frac{\hbar
v}{2a}\left[\left(\psi_{\boldsymbol{0}}^{1}\right)^{\dagger}\\!\\!\psi_{a\hat{x}}^{2}\\!-\\!i\left(\psi_{\boldsymbol{0}}^{1}\right)^{\dagger}\\!\\!\psi_{a\hat{y}}^{2}+\\!\left(\psi_{\boldsymbol{0}}^{1}\right)^{\dagger}\\!\\!\psi_{a\hat{z}}^{1}\\!\right.$
(27)
$\displaystyle\left.-\left(\psi_{\boldsymbol{0}}^{1}\right)^{\dagger}\\!\\!\psi_{-a\hat{x}}^{2}\\!\\!+\\!i\left(\psi_{\boldsymbol{0}}^{1}\right)^{\dagger}\\!\\!\psi_{-a\hat{y}}^{2}\\!\\!-\\!\left(\psi_{\boldsymbol{0}}^{1}\right)^{\dagger}\\!\\!\psi_{-a\hat{z}}^{1}\\!-\\!\text{h.c.}\right],$

while in the opposite case, it reads

$\displaystyle\mathcal{V}_{\\!2}\\!$ $\displaystyle=\\!-\frac{\hbar
v}{2a}\left[\left(\psi_{\boldsymbol{0}}^{2}\right)^{\dagger}\\!\\!\psi_{a\hat{x}}^{1}\\!+\\!i\left(\psi_{\boldsymbol{0}}^{2}\right)^{\dagger}\\!\\!\psi_{a\hat{y}}^{1}-\\!\left(\psi_{\boldsymbol{0}}^{2}\right)^{\dagger}\\!\\!\psi_{a\hat{z}}^{2}\\!\right.$
(28)
$\displaystyle\left.-\left(\psi_{\boldsymbol{0}}^{2}\right)^{\dagger}\\!\\!\psi_{-a\hat{x}}^{1}\\!\\!-\\!i\left(\psi_{\boldsymbol{0}}^{2}\right)^{\dagger}\\!\\!\psi_{-a\hat{y}}^{1}\\!+\\!\left(\psi_{\boldsymbol{0}}^{2}\right)^{\dagger}\\!\\!\psi_{-a\hat{z}}^{2}\\!-\\!\text{h.c.}\right].$

In either situation, the presence of bound-states and DoS deformations can be
investigated along the same lines of the full-vacancy, the only difference
being the projected perturbation matrix. More specifically, we write

$\mathcal{V}_{1/2}\\!=\\!\frac{\hbar
v}{2a}\left[\begin{array}[]{ccc}\mathbb{O}_{{\scriptscriptstyle
2\\!\times\\!2}}&\\!\\!\\!-\\!\boldsymbol{\Sigma}_{1/2}&\boldsymbol{\Sigma}_{1/2}\\\
\boldsymbol{\Sigma}_{1/2}^{T}&\mathbb{O}_{{\scriptscriptstyle
6\\!\times\\!6}}&\mathbb{O}_{{\scriptscriptstyle 6\\!\times\\!6}}\\\
\\!\\!\\!-\\!\boldsymbol{\Sigma}_{1/2}^{T}&\mathbb{O}_{{\scriptscriptstyle
6\\!\times\\!6}}&\mathbb{O}_{{\scriptscriptstyle
6\\!\times\\!6}}\end{array}\right],$ (29)

where the $(2\\!\times\\!6)$-dimensional $\boldsymbol{\Sigma}$–matrices are

$\displaystyle\boldsymbol{\Sigma}_{\text{u}}\\!$
$\displaystyle=\\!\left[\begin{array}[]{cccccc}0&1&0&-i&1&0\\\
0&0&0&0&0&0\end{array}\right]\,,$ (30c)
$\displaystyle\boldsymbol{\Sigma}_{\text{l}}\\!$
$\displaystyle=\\!\left[\begin{array}[]{cccccc}0&0&0&0&0&0\\\
1&0&i&0&0&-1\end{array}\\!\right].$ (30f)

The clean pGF is exactly the same as in Eq. (25), and therefore

$\det\left(\mathcal{I}\\!-\\!\overline{\boldsymbol{G}_{0}}(0)\\!\cdot\\!\mathcal{V}_{1/2}\right)\\!=\\!-\left(i\\!-\\!3g_{0}\right)^{2}$
(31)

which has a double root for $g_{0}\\!=\\!\nicefrac{{i}}{{3}}$. The
diagonalization of this operator confirms that its null-space is a two-
dimensional subspace generated by the removed orbital, plus a non-trivial
bound state that surrounds the vacant site. Similarly, the correction to the
DoS is exactly the same as Eq. (26) but with an added factor of
$\nicefrac{{1}}{{2}}$. The similar behavior between half- and full-vacancies
could have been anticipated by looking at a full-vacancy as a pair of half-
vacancies placed within the same unit cell. In the pGF formalism, these
correspond to local perturbations ($\mathcal{V}_{1}$ and $\mathcal{V}_{2}$)
that act in disjoint Hilbert subspaces not connected by the clean lattice
propagator. This prevents the two half-vacancies from hybridizing and their
resulting effects in the spectrum will be simply cumulative. For this reason,
we focus exclusively on full-vacancies in the remainder of this paper.

Figure 3: DoS correction due to a single (full) vacancy in an infinite lattice
(calculated with numerical differentiation of the lattice GF). The unperturbed
bandwidth of the model is marked by the outermost vertical dashed lines and,
unlike the conventional case, integrating $\delta\nu(E)$ over the entire band
yields $-2$ instead of zero.

##  III Microscopic Robustness of the Vacancy Bound States

After establishing the existence of vacancy-induced nodal bound states, we
move on to assess their robustness against additional disorder sources. For
that, we model the additional disordered landscape as an uncorrelated scalar
potential, $V_{d}(\boldsymbol{R})$. The Hamiltonian now reads

Figure 4: Probability density,
$\left|\Psi\\!\left(\mathbf{R}\right)\right|^{2}\\!\\!=\\!\left|\Psi_{1}\\!\left(\mathbf{R}\right)\right|^{2}\\!+\\!\left|\Psi_{2}\\!\left(\mathbf{R}\right)\right|^{2}$,
of the eigenstate closest to zero averaged over spherical shells of width $a$
centered on the vacancy. (a) Single central vacancy with no additional
disorder. Inset: The eigenvalues obtained near $E\\!=\\!0$ for finite samples
come in two flavors: those which scale to zero energy as
$L^{{\scriptscriptstyle-1}}$ (orange), akin to the predicted mean-level
spacing scaling around a clean Weyl node, and eigenvalues which scale faster
[as $\propto L^{{\scriptscriptstyle-3/2}}$] due to the more localized
character of the corresponding eigenstate (magenta). (b) Vacancy within three
random disordered environments.

$\mathcal{H}_{d}\\!=\\!\mathcal{H}_{0}\\!+\\!\mathcal{V}\\!+\\!\sum_{\mathbf{R}\in\mathcal{L}}V_{d}(\boldsymbol{R})\Psi_{\boldsymbol{R}}^{\dagger}\\!\cdot\\!\Psi_{\boldsymbol{R}}.$
(32)

To tackle this problem, we numerically diagonalize $\mathcal{H}_{d}$ around
$E\\!=\\!0$ using the scipy implementation of the implicitly restarted Lanczos
method [59, 60, 61]. Since the method converges better to non-clustered
eigenpairs in the borders of the spectrum, we apply it to
$\mathcal{H}_{d}^{2}$ instead and restrict the analysis to low-lying
eigenstates. Additionally, we consider cubic samples of side $L$, with a
single (full) vacancy at the center of each sample, supplemented by fixed
phase-twisted boundaries that open a finite-size gap
($\Delta_{f}\\!\propto\\!L^{{\scriptscriptstyle-1}}$) in the spectrum of
extended states. The nodal bound states will lie inside this finite-size gap,
as they are weakly affected by the boundary conditions.

In Fig. 4 a, we represent the radial distribution of the vacancy bound states
for different simulation sizes, in the absence of any additional disorder. The
results confirm the predicted zero-energy states with tails decaying as
$r^{{\scriptscriptstyle-2}}$. We note that the degeneracy of these states gets
slightly lifted by the boundary conditions but the corresponding eigenvalues
still tends to zero as $L^{{\scriptscriptstyle-\frac{3}{2}}}$, i.e., faster
than $\Delta_{f}$. Furthermore, the vacancy defect perturbs slightly the
extended states (now scattering states), which further expands the finite-size
gap 111This difference is too small to be seen in the inset of Fig. 4 and does
not affect the asymptotic scaling of $\Delta_{f}$ with the system size..

Next, we present an identical analysis with the vacancy’s surroundings endowed
with an uncorrelated random scalar potential uniformly drawn from
$[-\nicefrac{{W}}{{2}},\nicefrac{{W}}{{2}}]$. In principle, this alteration
dresses the LGF of the clean model, thus destroying most model-specific
symmetries. In Fig. 4 b, we present the radial wavefunctions of the two
eigenstates closest to zero energy obtained from the diagonalization of three
randomly generated disorder configurations. In all three cases, the states
feature the same normalizable power-law tail found in the clean case,
indicating that the vacancy bound states are indeed robust.

To further understand the effects of an uncorrelated disorder landscape, we
diagonalized $10^{{\scriptscriptstyle 4}}$ systems with randomly generated
disorder samples around a single central vacancy, focusing on determining the
four eigenpairs whose energies are the closest to the node. In addition to the
eigenenergies, we used the eigenfunctions to determine the inverse
participation ratio (IPR),

$\text{IPR}_{\Psi}\\!=\frac{\sum_{\mathbf{R}}\left(\left|\Psi_{\mathbf{R}}^{1}\right|^{2}\\!\\!+\\!\left|\Psi_{\mathbf{R}}^{2}\right|^{2}\right)^{2}}{\sum_{\mathbf{R}}\left|\Psi_{\mathbf{R}}^{1}\right|^{2}\\!\\!+\\!\left|\Psi_{\mathbf{R}}^{2}\right|^{2}},$
(33)

a simple quantity that allows one to distinguish well localized,
$\text{IPR}\\!\sim\\!\mathcal{O}(1)$, from delocalized states, for which
$\text{IPR}\\!\sim\\!\mathcal{O}(L^{{\scriptscriptstyle-3}})$. In Fig. 5 a and
5 b, we show histograms of the eigenenergies for three system sizes and two
disorder strengths, using twisted boundary conditions with a fixed twist angle
of $\nicefrac{{\pi}}{{3}}$ in all directions. These histograms borne out two
well-separated clusters formed by: (a) the two eigenstates closest to zero
energy which are broadened by disorder around $E\\!=\\!0$, but remain firmly
inside the finite-size gap, and (b) the ones corresponding to the largest
eigenvalues, being Bloch states that get perturbatively shifted towards the
node and broadened by disorder [28]. As confirmed by the corresponding IPRs,
the (a)-class states are strongly localized states which are still bound to
the central vacancy, while the (b)-class are disorder-dressed extended Weyl
states.

##  IV Quantum-Interference and Finite Concentration Effects

The previous results established that a single (full) vacancy defect gives
rise to a pair of zero-energy bound states with power-law-localized
wavefunctions. We now discuss the effect of coherent multiple scattering
events in realistic systems [51] which have a finite (nonzero) concentration
of point defects. The main question we ask here is whether the essential IPR
features of zero-energy states survive the unavoidable inter-vacancy
hybridization effects.

Our starting approach to this problem is based upon the exact diagonalization
of small systems. We consider WSM lattices with linear sizes up to $L=35$ and
a concentration ($n$) of randomly placed full-vacancies. By means of the
twisted boundary conditions, we open a finite-size gap that separates nodal
bound states from extended ones. The $2nL^{3}\\!+\\!4$ eigenpairs 222We expect
each vacancy to generate a pair of bound states, making this number
sufficiently large to capture most states within the finite-size gap. closest
to $E\\!=\\!0$ are then extracted using LD. In Fig. 5 c we represent a scatter
plot of the energies and corresponding IPRs of every eigenpair determined for
$2500$ random arrangements of vacancies with concentrations ranging from
$0.1\%$ to $1\%$ (per unit cell). The results clearly demonstrate that, in
spite of the proximity between vacancies, the system still features a large
number of high-IPR eigenstates which are flanked by a region of extended
states. Such a physical interpretation is clear from Fig. 5 d, where a 3D
bubble chart of $\left|\Psi(\mathbf{R})\right|^{2}$ is depicted for two
eigenstates randomly chosen from each of the regions.

Figure 5: (a) Histograms of the four energies closest to the Weyl node for a
vacancy surrounded by a random Anderson potential of strength
$W\\!=\\!0.01\hbar v/a$ (right) and $W\\!=\\!0.04\hbar v/a$ (left),
represented for systems with different linear sizes [$L=15$ (red), $L=25$
(blue) and $L=35$ (red)]. Black dashed lines indicate the clean energy levels
and, thus, the corresponding finite size gaps. (b) Scatter plot of the
energies and IPRs of the calculated eigenpairs. The two eigenstates closest to
zero energy are very localized in real-space, while the remaining pair is
clearly delocalized. (c) Energy-IPR scatter plot for the eigenstates closest
to the node obtained for $25000$ samples of a WSM with four concentrations of
randomly placed full-vacancies. Histograms of the IPRs are presented for each
total number of vacancies along the vertical axis. Different colors label
different vacancy numbers (here, the simulated cell size adjusted to guarantee
that $n$ is constant in all 4 cases presented). (d) Bubble chart
representation of the (squared) wavefunctions for two eigenstates picked at
random from the indicated regions: a heavily localized state around a few
vacancies, (left) and an extended state (right).

###  IV.1 Average Density of States

The LD study gives a qualitative picture of the structure of the eigenstates
surrounding a Weyl node, but its utility can be severely limited by loss of
spectral resolution, the finite number of eigenstates that are accessible and,
lastly, the attainable system sizes. Therefore, we now complement the LD
analysis with full-spectral simulations of the DoS of large systems by means
of the kernel polynomial method (KPM) [62]. As a first step, we present
results on the ensemble-averaged DoS for a large system with a linear size
$L\\!=\\!512$. This observable gives us information on how the spectrum is
modified by inter-vacancy hybridization effects, yielding a numerically exact
picture of the vacancy-induced resonances [27, 45] around the node. The KPM
calculations are carried out with domain decomposition and a stochastic trace
evaluation techniques as implemented in KITE [64]. The calculation employs
$M\\!=\\!2^{16}$ Chebyshev moments (corresponding to a spectral resolution
$\eta\\!=\\!10^{-4}\hbar v/a$), a Jackson kernel, and a sufficiently large
number of random vectors to yield highly converged results. Finally, the
results are averaged over random twisted boundary conditions which eliminates
the finite size mean-level spacing.

The average DoS obtained through the KPM is shown in Fig. 6. These high-
resolution results disclose a prominent spectral enhancement in and around the
node which indicates that, unlike what happens for ordinary on-site disorder
[42, 24, 23, 25, 46, 44, 26, 27, 28, 45], the DoS at $E\\!=\\!0$ gets quickly
lifted to a large value as $n$ increases. This pronounced effect is consistent
with the presence of robust nodal bound states, and validates the conclusions
of Sec. IV. Moreover, as the central DoS peak grows in height with increasing
$n$, a much wider symmetrical profile begins to emerge at its base. As shown
in Appendix B, the integral of this DoS correction is proportional to $n$,
which indicates that inter-vacancy hybridization is simply turning the bound
states of isolated (full) vacancies into scattering resonances within the
continuum. In Fig. 6 b, a closeup of this structure is shown, revealing a
finer comb-like structure of subsidiary peaks (sharp scattering resonances)
around the node for moderate defect concentrations ($n\lesssim 1\%$). These
peaks are more visible in Fig. 6 c, where their displacement as a function of
$n$ is also shown. The modulated structure in the DoS reported in this work is
a unique feature of 3D WSMs, which is absent in the analogous two-dimensional
problem [69, 63] (see inset to Fig. 6 a). The subsidiary peaks in the DoS are
robust to an additional weakly disordered potential, as is discussed in
Appendix B.

##  V Conclusion and Outlook

A combination of exact diagonalization and large-scale spectral methods allow
us to resolve the impact of point defects on the real-space electronic
structure of 3D $\mathcal{T}$-symmetric Weyl semimetals. Our results for a
lattice WSM model show that dilute concentrations of vacancies, a common
crystal imperfection, have a strong impact on the electronic properties in
stark contrast with uncorrelated on-site disorder models [42, 24, 31, 74, 23,
44, 45, 34, 46, 26, 27, 28, 25] which have been found to produce a minute
effect on the low-energy properties of Weyl systems. In fact, random vacancies
were shown to efficiently lift the nodal DoS thereby destabilizing the semi-
metallic phase even at very low concentrations. Moreover, we have also
demonstrated that quantum-interference effects between vacancies can yield a
peculiar modulated energy dependence of electronic observables which has no
analogue in two-dimensional Dirac systems [29, 70, 63].

Figure 6: DoS of a WSM lattice of linear size $L\\!=\\!512$ for selected
defect concentrations. (a) Bird’s-eye view of the DoS. The inset shows the
two-dimensional case for comparison. (b) Close-up of average DoS around
$\varepsilon\\!=\\!0$. (c) Subsidiary peaks in the DoS at low vacancy
concentration. Inset: Evolution of the the peak height ($P_{1},P_{2}\text{ and
}P_{3}$) with $n$.

While the average DoS displays a comb of subsidiary resonance peaks at finite
energies, we show in a companion paper [75] that the bulk dc-conductivity
mirrors this behavior through a series of sudden dips as the Fermi level is
varied. Therefore, upon tuning the carrier density in real samples [76] [or
even the defect concentration using H/He [53] or light-ion irradiation [54]],
we predict that bulk transport measurements will allow the observation of
interesting signatures of native point defects. These are expected to assume
chief importance in WSMs of the TaAs family, whose concentration of point
defects in high-quality crystals grown by chemical vapor transport is
experimentally known to be large [51]. At last, in Ref. [75] the authors
present a thorough investigation of physical consequences from these novel
effects of vacancy disorder in electronic structure of the Weyl nodes, with a
particular focus on experimentally accessible signatures from standard
transport and optical response measurements. The optical signatures have a
particular practical importance, as they do not rely on an external control
over the system’s Fermi energy.

###### Acknowledgements.

J.P.S.P., S.M.J, B.A. and J.M.V.P.L acknowledge support from the Portuguese
Foundation for Science and Technology (FCT) within the Strategic Funding
UIDB/04650/2020, and through projects No. POCI-01-0145-FEDER-028887 (J.P.S.P.,
S.M.J and J.M.V.P.L) and No. CEECIND/02936/2017 (B.A.). J.P.S.P. and S.M.J are
funded by FCT grants No. PD/BD/142774/2018 and PD/BD/142798/2018,
respectively. A.F. acknowledges support from the Royal Society (London)
through a Royal Society University Research Fellowship. The large-scale
calculations were undertaken on the HPC Viking Cluster of the University of
York. We thank J. Dieplinger, J.M.B Lopes dos Santos and A. Altland for
fruitful discussions. We also thank the anonymous referees for providing
valuable feedback.

##  Appendix A Lattice Green Function

Here, we outline the semi-analytical method employed to obtain the LGF with
arbitrary spectral resolution. As indicated in the main text, Eqs. (3)-(4),
the clean LGF is a $2\\!\times\\!2$ matrix that can be written as

$\boldsymbol{G}_{0}(\varepsilon,\\!\boldsymbol{\Delta
R})\\!=\\!\mathcal{I}_{\varepsilon}^{0}\\!\\!\left(\boldsymbol{\Delta
R}\right)\\!-\\!\\!\\!\\!\sum_{j=x,y,z}\\!\\!\\!\\!i\sigma^{j}\mathcal{I}_{\varepsilon}^{j}\\!\\!\left(\boldsymbol{\Delta
R}\right),$ (34)

with $\mathcal{I}_{\varepsilon}^{0,x,y,z}\left(n_{x},n_{y},n_{z}\right)$ being
four position- and energy-dependent integrals over the cubic fBZ. These
integrals are given by

$\displaystyle\mathcal{I}_{\varepsilon}^{0}\\!\left(\\!\boldsymbol{\Delta
R}\right)\\!$
$\displaystyle=\\!\\!\int_{\text{[-\leavevmode{p},\leavevmode{p}]}^{3}}\\!\\!\\!\\!\\!\\!\\!\\!\\!d^{{\scriptscriptstyle\text{(\\!3\\!)}}}\mathbf{k}\frac{\varepsilon\,e^{ik_{x}n_{x}}e^{ik_{y}n_{y}}e^{ik_{z}n_{z}}}{8\pi^{3}\left(\varepsilon^{2}+\left|\boldsymbol{\sin
k}\right|^{2}\right)}$ (35a)
$\displaystyle\mathcal{I}_{\varepsilon}^{x}\\!\left(\\!\boldsymbol{\Delta
R}\right)\\!$
$\displaystyle=\\!\\!\\!\int_{\text{[-\leavevmode{p},\leavevmode{p}]}^{3}}\\!\\!\\!\\!\\!\\!\\!\\!\\!d^{{\scriptscriptstyle\text{(\\!3\\!)}}}\mathbf{k}\frac{\sin\\!k_{x}e^{ik_{x}n_{x}}e^{ik_{y}n_{y}}e^{ik_{z}n_{z}}}{8\pi^{3}\left(\varepsilon^{2}+\left|\boldsymbol{\sin
k}\right|^{2}\right)},$ (35b)

For a single-vacancy calculation in the lattice WSM, the finite support of the
perturbation dictates that only three of these integrals are required, namely:

$\displaystyle\mathcal{I}_{\varepsilon}^{0}\\!\left(0,\\!0,\\!0\right)\\!$
$\displaystyle=\\!\int_{\text{[-\leavevmode{p},\leavevmode{p}]}^{3}}\\!\\!\\!\\!\\!\\!\\!\\!\\!d^{{\scriptscriptstyle\text{(\\!3\\!)}}}\mathbf{k}\frac{\varepsilon}{8\pi^{3}\left(\varepsilon^{2}\\!+\\!\left|\boldsymbol{\sin
k}\right|^{2}\right)}$ (36a)
$\displaystyle\mathcal{I}_{\varepsilon}^{x}\\!\left(1,\\!0,\\!0\right)\\!$
$\displaystyle=\\!\int_{\text{[-\leavevmode{p},\leavevmode{p}]}^{3}}\\!\\!\\!\\!\\!\\!\\!\\!\\!d^{{\scriptscriptstyle\text{(\\!3\\!)}}}\mathbf{k}\frac{\sin\\!k_{x}e^{ik_{x}}}{8\pi^{3}\left(\varepsilon^{2}\\!+\\!\left|\boldsymbol{\sin
k}\right|^{2}\right)}$ (36b)
$\displaystyle\mathcal{I}_{\varepsilon}^{0}\\!\left(2,\\!0,\\!0\right)\\!$
$\displaystyle=\\!\int_{\text{[-\leavevmode{p},\leavevmode{p}]}^{3}}\\!\\!\\!\\!\\!\\!\\!\\!\\!d^{{\scriptscriptstyle\text{(\\!3\\!)}}}\mathbf{k}\frac{\varepsilon
e^{2ik_{x}}}{8\pi^{3}\left(\varepsilon^{2}\\!+\\!\left|\boldsymbol{\sin
k}\right|^{2}\right)},$ (36c)

which define three complex-valued functions $f_{\varepsilon}$,
$g_{\varepsilon}$ and $h_{\varepsilon}$, respectively. In all previous
integrals $\varepsilon$ is to be taken as a complex number with a positive
infinitesimal imaginary part (guaranteeing that the LGFs are retarded). In all
three cases, analytical progress can be made by first considering the one-
dimensional integral,

$I_{1}\\!\left(z\right)\\!=\\!\int_{-\pi}^{\pi}\\!\\!\\!\\!\\!du\frac{1}{2\pi\left(z\\!+\\!\sin^{2}\\!u\right)}$
(37)

where $z\in\mathbb{C}$. This integral can be solved by standard contour
integration in the complex variable $w\\!=\\!\exp\\!\left(iu\right)$, yielding

$I_{1}\\!(z)\\!=\\!\frac{1}{\sqrt{z}\sqrt{z-1}}.$ (38)

Therefore,

$I_{1}^{\pm}\\!\\!\left(x\right)\\!=\\!\\!\\!\lim_{\eta\to
0^{{\scriptscriptstyle\pm}}}\\!\left[I_{1}\\!\\!\left(x\\!+\\!i\eta\right)\right]\\!=\\!\begin{cases}\frac{\text{sgn}(x)}{\sqrt{x(x\\!-\\!1)}}&\\!\\!\\!x\\!\notin\\!\left[0,1\right]\\\
\mp\frac{i}{\sqrt{x(1\\!-\\!x)}}&\\!\\!\\!x\\!\in\\!\left[0,1\right]\end{cases}.$
(39)

From Eq. (39), it is easy to recognize that all three integrals in Eqs.
(36a)-(36c) can be written as two-dimensional $\mathbf{k}$-integrals involving
$I_{1}^{{\scriptscriptstyle\pm}}\\!(x)$. Hence, they read

Figure 7: Numerical integrals that define the complex-valued functions
$f_{\varepsilon}$, $g_{\varepsilon}$ and $h_{\varepsilon}$. The node is placed
at $\varepsilon\\!=\\!0$. Van Hove singularities (derivative discontinuities)
are marked as vertical dimmed lines. All the curves were obtained in an
regular energy mesh of spacing $10^{{\scriptscriptstyle-5}}\hbar v/a$.

$\displaystyle\\!\\!f_{E\pm i0^{+}}\\!$
$\displaystyle=\\!\frac{E}{4\pi^{2}}\int_{\text{[-\leavevmode{p},\leavevmode{p}]}^{2}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!d^{{\scriptscriptstyle\text{(\\!2\\!)}}}\mathbf{k}\,I_{1}^{\pm}\\!\left(E^{2}\\!\\!+\\!\sin^{2}\\!k_{x}\\!+\\!\sin^{2}\\!k_{y}\right)$
(40a) $\displaystyle\\!\\!g_{E\pm i0^{+}}\\!$
$\displaystyle=\\!\frac{1}{4\pi^{2}}\\!\\!\int_{\text{[-\leavevmode{p},\leavevmode{p}]}^{2}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!d^{{\scriptscriptstyle\text{(\\!2\\!)}}}\mathbf{k}\sin\\!k_{x}e^{ik_{x}}\,I_{1}^{\pm}\\!\left(E^{2}\\!\\!+\\!\sin^{2}\\!k_{x}\\!+\\!\sin^{2}\\!k_{y}\right)$
(40b) $\displaystyle\\!\\!h_{E\pm i0^{+}}\\!$
$\displaystyle=\\!\frac{E}{4\pi^{2}}\int_{\text{[0,\leavevmode{p}]}^{2}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!d^{{\scriptscriptstyle\text{(\\!2\\!)}}}\mathbf{k}e^{2ik_{x}}\,I_{1}^{\pm}\\!\left(E^{2}\\!\\!+\\!\sin^{2}\\!k_{x}\\!+\\!\sin^{2}\\!k_{y}\right),$
(40c)

where now $E$ is the (real-valued) energy parameter and the $\pm$ denotes the
sign of the energy. These two-dimensional integrals can be numerically
evaluated for an arbitrarily fine mesh of $E$. As $\eta$ was formally taken to
zero, the spectral resolution is only limited by the spacing of this mesh.
Moreover, using the symmetries — $f_{-E\pm i0^{+}}\\!=\\!-f_{E\pm
i0^{+}}^{*}$, $g_{-E\pm i0^{+}}\\!=\\!g_{E\pm i0^{+}}^{*}$, $h_{-E\pm
i0^{+}}\\!=\\!-h_{E\pm i0^{+}}^{*}$ — it is enough to evaluate the said
integrals for $E\\!>\\!0$. The results are shown in Fig. 7.

##  Appendix B Additional Numerical Results

Here, we complement the KPM results shown in Sec.IV.1 for the mean DoS in the
presence of a finite concentration of vacancies. These results serve as a
support for some of the claims made in the main text.

First, we analyze the region of enhanced mean DoS around the nodal energy
($E\\!=\\!0$). For a very low vacancy concentration, the sole feature is a
sharp peak at the node, which amounts to a cumulative contribution of all
single-vacancy bound states to the intensive DoS. As the concentration
increases, one starts to see more structure, in the form of a broadened base
of this peak (as highlighted in the main panel of Fig. 8 a). We attribute this
to a progressive hybridization of states bound to nearby vacancies that lifts
the degeneracies away from the node, thus turning these states into long-
living resonances. To corroborate this idea, we analyze the dependence of the
integrated change in the DoS within this central region. Broadly speaking,
this quantity represents the number of states introduced near the node by the
vacancies, per unit volume. In the inset to Fig. 8 a, we show that the
integrated DoS change induced by the vacancies scales exactly as
$n/(100\\!-\\!n)$. This is consistent with the picture in which each missing
Wannier state introduces exactly one eigenstate in the node and one eigenstate
around the node.

Another point concerns with the robustness of the features of DoS to
microscopic details of the underlying lattice model. As we have done for a
single (full) vacancy, here also we probe this robustness by introducing an
additional Anderson potential (of strength $W$). The KPM results for the mean
DoS are shown in Fig. 8 b. Clearly, the main features of the DoS, i.e., the
central enhancement and the subsidiary peaks, remain untouched for suitably
small $W$. This supports the claim that our results will hold for a wide range
of WSM systems.

Figure 8: (a) Full integral of the correction do the mean DoS as a function of
the vacancy concentration. Inset: Integral of
$\delta\nu(\varepsilon)\\!=\\!\nu(\varepsilon)\\!-\\!\nu_{0}(\varepsilon)$
within the region represented plotted in the main panel. (b) Mean DoS
calculated for two concentrations of vacancies and an additional potential
disorder.

## References

  * Fu and Kane [2007] L. Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Murakami [2007] S. Murakami, New J. Phys. 9, 356 (2007).
  * Halász and Balents [2012] G. B. Halász and L. Balents, Phys. Rev. B 85, 035103 (2012).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * Weyl [1929] H. Weyl, Proc. Natl. Acad. Sci. U.S.A. 15, 323 (1929).
  * Wan  _et al._ [2011] X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Phys. Rev. B 83, 205101 (2011).
  * Xu  _et al._ [2011] G. Xu, H. Weng, Z. Wang, X. Dai, and Z. Fang, Phys. Rev. Lett. 107, 186806 (2011).
  * Burkov and Balents [2011] A. A. Burkov and L. Balents, Phys. Rev. Lett. 107, 127205 (2011).
  * Koshino and Hizbullah [2016] M. Koshino and I. F. Hizbullah, Phys. Rev. B 93, 045201 (2016).
  * Berry [1984] M. V. Berry, Proc. R. Soc. Lond. A. Math. Phys. Sci. , 45 (1984).
  * Hosur [2012] P. Hosur, Phys. Rev. B 86, 195102 (2012).
  * Witten [2016] E. Witten, Riv. del Nuovo Cimento 39, 313 (2016).
  * [14] F. D. M. Haldane, arXiv:1401.0529 .
  * Hashimoto  _et al._ [2017] K. Hashimoto, T. Kimura, and X. Wu, Prog. Theor. Exp. Phys. 2017, 053I01 (2017).
  * Adler [1969] S. L. Adler, Phys. Rev. 177, 2426 (1969).
  * Bell and Jackiw [1969] J. S. Bell and R. Jackiw, Il Nuovo Cimento A (1965-1970) 60, 47 (1969).
  * Son and Spivak [2013] D. T. Son and B. Z. Spivak, Phys. Rev. B 88, 104412 (2013).
  * Zhang  _et al._ [2016] C.-L. Zhang, S.-Y. Xu, I. Belopolski, Z. Yuan, Z. Lin, B. Tong, G. Bian, N. Alidoust, C.-C. Lee, S.-M. Huang, T.-R. Chang, G. Chang, C.-H. Hsu, H.-T. Jeng, M. Neupane, D. Sanchez, H. Zheng, J. Wang, H. Lin, C. Zhang, H.-Z. Lu, S.-Q. Shen, T. Neupert, M. Zahid Hasan, and S. Jia, Nat. Comm. 7, 10735 (2016).
  * Burkov [2017] A. A. Burkov, Phys. Rev. B 96, 041110(R) (2017).
  * Nandy  _et al._ [2017] S. Nandy, G. Sharma, A. Taraphder, and S. Tewari, Phys. Rev. Lett. 119, 176804 (2017).
  * Nielsen and Ninomiya [1983] H. Nielsen and M. Ninomiya, Physics Letters B 130, 389 (1983).
  * Bera  _et al._ [2016] S. Bera, J. D. Sau, and B. Roy, Phys. Rev. B 93, 201302(R) (2016).
  * Pixley  _et al._ [2015] J. H. Pixley, P. Goswami, and S. Das Sarma, Phys. Rev. Lett. 115, 076601 (2015).
  * Pixley  _et al._ [2016] J. H. Pixley, P. Goswami, and S. Das Sarma, Phys. Rev. B 93, 085103 (2016).
  * Wilson  _et al._ [2020] J. H. Wilson, D. A. Huse, S. Das Sarma, and J. H. Pixley, Phys. Rev. B 102, 100201(R) (2020).
  * Santos Pires  _et al._ [2021] J. P. Santos Pires, B. Amorim, A. Ferreira, İ. Adagideli, E. R. Mucciolo, and J. M. Viana Parente Lopes, Phys. Rev. Research 3, 013183 (2021).
  * Pixley and Wilson [2021] J. Pixley and J. H. Wilson, Ann. Phys. , 168455 (2021).
  * Fradkin [1986] E. Fradkin, Phys. Rev. B 33, 3263 (1986).
  * Roy and Das Sarma [2014] B. Roy and S. Das Sarma, Phys. Rev. B 90, 241112(R) (2014).
  * Syzranov  _et al._ [2015a] S. V. Syzranov, L. Radzihovsky, and V. Gurarie, Phys. Rev. Lett. 114, 166601 (2015a).
  * Slager  _et al._ [2017] R.-J. Slager, V. Juričić, and B. Roy, Phys. Rev. B 96, 201401(R) (2017).
  * Shen  _et al._ [2020] J. Shen, Q. Yao, Q. Zeng, H. Sun, X. Xi, G. Wu, W. Wang, B. Shen, Q. Liu, and E. Liu, Phys. Rev. Lett. 125, 086602 (2020).
  * Syzranov and Radzihovsky [2018] S. V. Syzranov and L. Radzihovsky, Ann. Rev. Cond. Mat. Phys. 9, 35 (2018).
  * Roy  _et al._ [2018] B. Roy, R.-J. Slager, and V. Juričić, Phys. Rev. X 8, 031076 (2018).
  * Janssen [1994] M. Janssen, Int. J. Mod. Phys. B 8, 943 (1994).
  * Evers and Mirlin [2008] F. Evers and A. D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008).
  * Shindou and Murakami [2009] R. Shindou and S. Murakami, Phys. Rev. B 79, 045321 (2009).
  * Goswami and Chakravarty [2011] P. Goswami and S. Chakravarty, Phys. Rev. Lett. 107, 196803 (2011).
  * Syzranov  _et al._ [2015b] S. V. Syzranov, V. Gurarie, and L. Radzihovsky, Phys. Rev. B 91, 035133 (2015b).
  * Kobayashi  _et al._ [2014] K. Kobayashi, T. Ohtsuki, K.-I. Imura, and I. F. Herbut, Phys. Rev. Lett. 112, 016402 (2014).
  * Nandkishore  _et al._ [2014] R. Nandkishore, D. A. Huse, and S. L. Sondhi, Phys. Rev. B 89, 245110 (2014).
  * Gurarie [2017] V. Gurarie, Phys. Rev. B 96, 014205 (2017).
  * Buchhold  _et al._ [2018a] M. Buchhold, S. Diehl, and A. Altland, Phys. Rev. Lett. 121, 215301 (2018a).
  * Buchhold  _et al._ [2018b] M. Buchhold, S. Diehl, and A. Altland, Phys. Rev. B 98, 205134 (2018b).
  * Ziegler and Sinner [2018] K. Ziegler and A. Sinner, Phys. Rev. Lett. 121, 166401 (2018).
  * Q. Lv  _et al._ [2021] B. Q. Lv, T. Qian, and H. Ding, Rev. Mod. Phys. 93, 025002 (2021).
  * Ghimire  _et al._ [2015] N. J. Ghimire, Y. Luo, M. Neupane, D. J. Williams, E. D. Bauer, and F. Ronning, J. Phys.: Condens. Matt. 27, 152201 (2015).
  * Buckeridge  _et al._ [2016] J. Buckeridge, D. Jevdokimovs, C. R. A. Catlow, and A. A. Sokol, Phys. Rev. B 94, 180101(R) (2016).
  * Yu and Zhang [2016] Z. G. Yu and Y.-W. Zhang, Phys. Rev. B 94, 195206 (2016).
  * Besara  _et al._ [2016] T. Besara, D. A. Rhodes, K.-W. Chen, S. Das, Q. R. Zhang, J. Sun, B. Zeng, Y. Xin, L. Balicas, R. E. Baumbach, E. Manousakis, D. J. Singh, and T. Siegrist, Phys. Rev. B 93, 245152 (2016).
  * Liu  _et al._ [2016] H. W. Liu, P. Richard, L. X. Zhao, G.-F. Chen, and H. Ding, J. Phys.: Condens. Matt. 28, 295401 (2016).
  * Zhang  _et al._ [2019] Z.-J. Zhang, Y.-L. Fu, W. Cheng, and F.-S. Zhang, Computational Materials Science 160, 9 (2019).
  * Fu  _et al._ [2020] Y.-L. Fu, H.-B. Sang, W. Cheng, and F.-S. Zhang, Mater. Today Commun. 24, 100939 (2020).
  * Lehtinen  _et al._ [2010] O. Lehtinen, J. Kotakoski, A. V. Krasheninnikov, A. Tolvanen, K. Nordlund, and J. Keinonen, Phys. Rev. B 81, 153401 (2010).
  * Lin  _et al._ [2016] Z. Lin, B. R. Carvalho, E. Kahn, R. Lv, R. Rao, H. Terrones, M. A. Pimenta, and M. Terrones, 2D Mater. 3, 022002 (2016).
  * Bertoldo  _et al._ [2022] F. Bertoldo, S. Ali, S. Manti, and K. S. Thygesen, npj Computational Materials 8, 56 (2022).
  * Xing  _et al._ [2020] Y. Xing, J. Shen, H. Chen, L. Huang, Y. Gao, Q. Zheng, Y.-Y. Zhang, G. Li, B. Hu, G. Qian, L. Cao, X. Zhang, P. Fan, R. Ma, Q. Wang, Q. Yin, H. Lei, W. Ji, S. Du, H. Yang, W. Wang, C. Shen, X. Lin, E. Liu, B. Shen, Z. Wang, and H.-J. Gao, Nat Commun 11, 5613 (2020).
  * Lanczos [1950] C. Lanczos, J. Res. Natl. Inst. Stand. Technol. 45, 255 (1950).
  * Paige [1980] C. Paige, Linear Algebra Appl. 34, 235 (1980).
  * Lehoucq and Sorensen [1996] R. Lehoucq and D. C. Sorensen, SIAM J. Matrix Anal. Appl 17, 789 (1996).
  * Weiße  _et al._ [2006] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, Rev. Mod. Phys. 78, 275 (2006).
  * Ferreira and Mucciolo [2015] A. Ferreira and E. R. Mucciolo, Phys. Rev. Lett. 115, 106601 (2015).
  * João  _et al._ [2020] S. M. João, M. Anđelković, L. Covaci, T. G. Rappoport, J. M. Viana Parente Lopes, and A. Ferreira, R. Soc. open Sci. 7, 191809 (2020).
  * Siegel [1978] R. Siegel, Journal of Nuclear Materials 69-70, 117 (1978).
  * Crawford and Slifkin [1972] J. H. Crawford and L. M. Slifkin,  _Point Defects in Solids_ (Springer, Berlin, Germany, 1972).
  * Schaefer  _et al._ [1999] H.-E. Schaefer, K. Frenner, and R. Würschum, Phys. Rev. Lett. 82, 948 (1999).
  * Mahan [1995] G. D. Mahan, Int. J. Mod. Phys. B 9, 1313 (1995).
  * Pereira  _et al._ [2006] V. M. Pereira, F. Guinea, J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 96, 036801 (2006).
  * Ostrovsky  _et al._ [2014] P. M. Ostrovsky, I. V. Protopopov, E. J. König, I. V. Gornyi, A. D. Mirlin, and M. A. Skvortsov, Phys. Rev. Lett. 113, 186803 (2014).
  * Chiu  _et al._ [2016] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * Note [1] This difference is too small to be seen in the inset of Fig.+.1667em4 and does not affect the asymptotic scaling of $\Delta_{f}$ with the system size.
  * Note [2] We expect each vacancy to generate a pair of bound states, making this number sufficiently large to capture most states within the finite-size gap.
  * Syzranov  _et al._ [2016] S. Syzranov, V. Gurarie, and L. Radzihovsky, Ann. Phys. 373, 694 (2016).
  * Santos Pires  _et al._ [2022] J. P. Santos Pires, S. M. João, A. Ferreira, B. Amorim, and J. M. Viana Parente Lopes, arXiv:cond-mat/2205.15123 (2022).
  * Nishihaya  _et al._ [2018] S. Nishihaya, M. Uchida, Y. Nakazawa, M. Kriener, Y. Kozuka, Y. Taguchi, and M. Kawasaki, Sci. Adv. 4, eaar5668 (2018).