note-name = , use-sort-key = false

# Supplemental material

##  I Basic considerations

The larger the supercell gets, the more bands it contains in the given
frequency range. This effect is known as band-folding and it is a fundamental
property of the Brillouin zone. Several ways of "unfolding" the band structure
have been proposed in the literature, see, e.g., Refs. [1, 2]. Nevertheless,
it is still very hard to draw direct relations between bands from different
supercells as the band-folding effects seem to also depend on the confinement
dimensionality of the bands. Our approach circumvents this problem by not
requiring such direct relations between the folded and unfolded bands. Namely,
in the $\widetilde{V}^{\alpha}/\widetilde{E}$ plots for large enough
$N>N_{0}$, each band of the investigated $N$ supercell is clearly either above
or below all the bands with similar frequencies of the reference $N_{0}$
supercell.

Our method can be viewed as an extension of the multifractality techniques,
which are usually formulated for scaling of the participation ratio [3, 4, 5].
In the first step of our framework we also perform a similar scaling for the
mode volume in Eq. (LABEL:main-eq:vmodeecavD). However, by limiting oneself to
only this scaling behavior, it would require very large supercell size $N$ to
draw conclusive results for the whole spectrum, as described below and in Fig.
S-1. To complicate the matter even more, due to the band folding effects,
translating the results obtained for this large $N$ to experimentally
interesting supercells of moderate size would be extremely complicated, if at
all possible. By subsequently introducing the confinement energy as another
analyzed quantity with different scaling exponent than the mode volume and the
auxiliary power parameter $\alpha$, our method allows for determination of
confinement throughout the whole spectrum for relatively low supercell sizes
and, as described in the previous paragraph, free of complications related to
band folding.

In our definition of $E_{\mathrm{C}}$ in (LABEL:main-ecavdef) we choose the
volume $V_{\mathrm{C}}$ to be centered around the cavity, where the
confinement is expected. Nevertheless, this is not a strict requirement.
Choosing the volume $V_{\mathrm{C}}$ at a different location within the
supercell will still yield correct scaling results, albeit possibly affecting
the convergence of the method and thus requiring larger supercell sizes
$N,N_{0}$ to properly classify the whole spectrum. It is however necessary,
after $V_{\mathrm{C}}$ is selected, to keep this volume constant throughout
the scaling, so that the exponent $c-D$ appears in the scaling relation for
$\widetilde{E}$ in (LABEL:main-eq:vmodeecavD), since this is crucial for being
able to tune the values of the critical exponent $\kappa$ as needed.

The dimension of the supercell used to model the system may differ from the
system dimensionality $D$ in some cases. For example, in 2D electronic
materials $D=2$, but the charge density usually extends slightly into the
third dimension and, moreover, can be zero in the 2D plane where the atomic
nuclei are positioned. This case is therefore traditionally modeled by a 3D
supercell and thus the integration in Eqs. (LABEL:main-mvdef) and (LABEL:main-
ecavdef) must be performed over its whole three-dimensional volume, i.e., the
volumes $V_{\mathrm{S}}$ and $V_{\mathrm{C}}$ will be three-dimensional.
Nevertheless, since the structure of the material is still two-dimensional,
the subsequent scaling analysis of confinement can be performed in $D=2$
dimensions with the third dimension kept constant, i.e., adding unit cells
only along the planar material, as we did in the example of a hexagonal BN
with N vacancy.

In this Letter we are only concerned with integer defect dimensionalities $d$
in systems of dimension $D\leq 3$, as these are the most relevant for
experiments. It is however straightforward to extend our technique to also
include fractal defects with fractional dimensionalities and higher
dimensional spaces. One can directly apply the derivations in this
Supplemental material to find the optimal values of the auxiliary power
$\alpha$ in case of fractional $d$ or arbitrarily high $D$. These
possibilities further emphasize the generality of our scaling theory.

##  II Sparsity measures: Analogies of mode volume and participation ratio

In our confinement analysis, we employ the mode volume $V_{\mathrm{M}}$ as one
of two key quantities, defined in Eq. (LABEL:main-mvdef). This is, however,
not the only possible choice, and a whole class of quantities can be
substituted for $V_{\mathrm{M}}$ in our confinement identification method.

Ref. [6] is the first to rigorously define the so-called measures of sparsity.
In the most simple terms, measures of sparsity are functions quantifying the
inequality of distribution of values within a given set. One such measure, for
the set of values $\\{m_{1},\ldots,m_{M}\\}\in\mathbb{R}^{M}$, is a $pq$-mean:

$\mu_{p,q}\colonequals\frac{\left(\frac{1}{M}\sum\limits_{m=1}^{M}c_{j}^{p}\right)^{\nicefrac{{1}}{{p}}}}{\left(\frac{1}{M}\sum\limits_{m=1}^{M}c_{j}^{q}\right)^{\nicefrac{{1}}{{q}}}}.$
(S.1)

Ref. [6] has shown that, for $p<q<\infty$, $\mu_{p,q}$ is indeed a viable
measure of sparsity.

There is clear correspondence between the inequality of $W(\mathbf{x})$ over
the volume $V_{\mathrm{S}}$ for a given band and its confinement: For a
confined band, $W(\mathbf{x})$ is very high in certain parts of
$V_{\mathrm{S}}$ but low everywhere else, while for an extended band
$W(\mathbf{x})$ is approximately equally distributed over the whole
$V_{\mathrm{S}}$. To see this quantitatively, we write the straightforward
continuous generalization of (S.1) to quantify the inequality of
$W(\mathbf{x})$ over the volume $V_{\mathrm{S}}$ as

$U_{p,q}\colonequals\frac{\left(\frac{1}{V_{\mathrm{S}}}\int\limits_{V_{\mathrm{S}}}W^{p}(\mathbf{x})\mathrm{d}{V}\right)^{\nicefrac{{1}}{{p}}}}{\left(\frac{1}{V_{\mathrm{S}}}\int\limits_{V_{\mathrm{S}}}W^{q}(\mathbf{x})\mathrm{d}{V}\right)^{\nicefrac{{1}}{{q}}}}.$
(S.2)

Our normalized mode volume $\widetilde{V}$, as defined by Eq. (LABEL:main-
venorm), turns out to be a continuous version of the $pq$-mean with $p=1$ and
$q=\infty$, i.e., $\widetilde{V}=U_{1,\infty}$. The fact that $q=\infty$ means
that $\widetilde{V}$ has slightly weaker properties than the sparsity measure
defined by Ref. [6]. This detail is, however, unimportant for our confinement
identification purposes.

For the purpose of confinement analysis, one can substitute any other
$pq$-mean $U_{p,q}$ for $\widetilde{V}$ and our method will still work
properly, albeit possibly with different scaling power than we obtained for
$\widetilde{V}$ in Eq. (LABEL:main-eq:vmodeecavD). Specifically, substituting
$p=1$, $q=2$ yields

$U_{1,2}\colonequals\frac{\int\limits_{V_{\mathrm{S}}}W(\mathbf{x})\mathrm{d}{V}}{\sqrt{V_{\mathrm{S}}}\sqrt{\int\limits_{V_{\mathrm{S}}}W^{2}(\mathbf{x})\mathrm{d}{V}}}.$
(S.3)

The quantity $U_{1,2}$ corresponds to the square root of normalized
participation ratio:

$U_{1,2}=\sqrt{\frac{P}{V_{\mathrm{S}}}},$ (S.4)

where the participation ratio $P$ [7, 8, 9] is given by

$P\colonequals\frac{\left(\int_{V_{\mathrm{S}}}W(\mathbf{x})\mathrm{d}{V}\right)^{2}}{\int_{V_{\mathrm{S}}}W^{2}(\mathbf{x})\mathrm{d}{V}}.$
(S.5)

Finally, since the $pq$-mean is a non-negative dimensionless quantity, its
various powers are equivalent in characterizing confinement and thus one can
use $U_{1,2}^{2}$ instead of $U_{1,2}$. This illustrates that, in our
approach, instead of normalized mode volume $\widetilde{V}$, normalized
participation ratio can be used as well, similarly to any other $pq$-mean
$U_{p,q}$.

##  III Failure of multifractality

In the main text we claim that our method surpasses multifractality analysis
because it is able to identify bands in much smaller supercells than is
possible by multifractality. We illustrate this by applying the simple
multifractality scaling to our example of 3D inverse woodpile photonic crystal
with two linear defects used in the main text. We analyze the $N=4$ supercell
and use the $N_{0}=2$ supercell as a reference. We scale either the mode
volume $\widetilde{V}$, in Fig. S-1(a), or the confinement energy
$\widetilde{E}$, in Fig. S-1(b).

From Fig. S-1(a), we observe that for some bands the value of $\widetilde{V}$
decreases significantly as the supercell grows and for some it stays roughly
constant. One would immediately expect the bands with constant $\widetilde{V}$
to be extended, according to Eq. (LABEL:main-eq:vmodeecavD). However,
differentiating between the $c=2$ and $c=3$ bands is basically impossible.
Moreover, even the distinction between the "decreasing" and "constant" mode
volume can be rather unclear, for example around the frequencies
$\tilde{\omega}\approx 0.5$.

Similarly, in Fig. S-1(b), for some bands the value of $\widetilde{E}$
decreases significantly as the supercell grows and for some it stays roughly
constant. In this case, the distinction between the various confinement
dimensionalities is even blurrier than for the mode volume.

As it is suggestively implied by the combined Fig. S-1, looking at the scaling
of both $\widetilde{V}$ and $\widetilde{E}$ at the same time and comparing
these yields more information than analyzing only one of them as does the
traditional multifractality approach. This already brings us very close to our
confinement analysis method as described in the main text, which represents a
systematic approach to the simultaneous scaling of both $\widetilde{V}$ and
$\widetilde{E}$.

Figure S-1:  Failure of the standard multifractality approach to identify
confined bands for small supercell sizes. Each band is represented by a point.
(a) Scaling of the normalized mode volume. (b) Scaling of the normalized
confinement energy. For some bands the plotted values change significantly
when scaling from $N_{0}=2$ to $N=4$ while for others they remain
approximately constant. This variety of behaviors indicates that different
confinement dimensionalities $c$ are present in the band spectrum, as
predicted by the multifractality. Nevertheless, from each of these plots on
its own, it is not possible to clearly assign confinement dimensionalities to
specific bands and there is thus a need to move beyond the traditional
multifractality method. One can, albeit, infer interesting information on
confinement properties upon analyzing both plots simultaneously, which already
brings us close to the essence of our technique of scaling the ratio of
$\widetilde{V}^{\alpha}/\widetilde{E}$.

##  IV Derivation of the scaling equation

In a homogeneous $D$-dimensional system, waves propagate freely in all
directions. Adding a periodic structure with no defects may restrict some ways
of propagation, but the waves will still be extended, i.e., they can propagate
from any given unit cell to any other one. Upon introducing a defect of
dimensionality $d\leq D$, some waves may couple to this defect and only
propagate within it, thus never being able to achieve the non-defect unit
cells. A wave coupled to a defect of dimension $d$ has confinement
dimensionality $c=D-d$, i.e., it is confined in $c$ dimensions. This effect is
very strongly observable upon scaling of the system size. To analyze the
confinement dimensionality of waves, one can utilize the scaling of mode
volume and confinement energy.

Mathematically, the behavior of mode volume, defined by (LABEL:main-mvdef),
for such a wave is given by

$\displaystyle V_{\mathrm{M}}$ $\displaystyle=AN^{d}+\mathcal{O}(N^{d-1})$
(S.6) $\displaystyle=AN^{D-c}+\mathcal{O}(N^{D-c-1}),$

where $A$ is a constant independent of $N$. In fact, $A$ is actually a
functional depending on the permittivity distribution $\epsilon(\mathbf{x})$,
the wave-vector $\mathbf{k}$ and the frequency $\omega$. Eq. (S.6) holds
beyond the localization length of a given wave and expresses the fact that the
mode volume of a confined wave can only grow within the geometrical
constraints of the defect, with a small contribution of decaying waves in
other directions. Upon normalization by $V_{\mathrm{S}}=N^{D}$, we obtain

$\widetilde{V}=AN^{-c}+\mathcal{O}(N^{-c-1}).$ (S.7)

Similarly, the normalized confinement energy, defined by (LABEL:main-ecavdef)
and (LABEL:main-eq:vmodeecavD), behaves as

$\displaystyle\widetilde{E}$ $\displaystyle=BN^{-d}+\mathcal{O}(N^{-d-1})$
(S.8) $\displaystyle=BN^{c-D}+\mathcal{O}(N^{c-D-1}),$

where $B$ is a constant independent of $N$. Analogously to $A$, $B$ is in fact
a functional depending on the permittivity distribution
$\epsilon(\mathbf{x})$, the wave-vector $\mathbf{k}$ and the frequency
$\omega$. Eq. (S.8) holds beyond the localization length of the given wave and
expresses the fact that the confinement energy of a confined wave can only
escape from the volume $V_{0}$ in ways that obey the geometrical constraints
of the defect. There can also be additional, fast decaying, leakage expressed
by the term $\mathcal{O}(N^{c-D-1})$.

We combine the Eqs. (S.7) and (S.8) and, upon neglecting the big-$\mathcal{O}$
contributions, we obtain for a wave with the confinement dimensionality $c$ in
the limit $N\rightarrow\infty$:

$\frac{\widetilde{V}^{\alpha}}{\widetilde{E}}=\frac{\left(AN^{-c}\right)^{\alpha}}{BN^{c-D}}=CN^{-(\alpha+1)c+D},$
(S.9)

where $C=A^{\alpha}/B$. Eq. (S.9) represents exactly our scaling relation in
Eq. (LABEL:main-eq:scaling).

##  V Determining the auxiliary powers $\alpha$

In our method, we are able to distinguish between the bands with $c<j$ and
$c\geq j$ by tracking the behavior of $\widetilde{V}^{\alpha}/\widetilde{E}$
according to Eq. (LABEL:main-eq:scaling), as the supercell size changes from a
smaller reference size $N_{0}$ to the investigated size $N$. By varying $j$
over all the values $0<j\leq D$, this approach allows us to completely
classify all wave bands in the spectrum based on their confinement
dimensionality. For this approach to work, the power $\alpha$ has to be chosen
so that $\kappa<0$ for all $c\geq j$ and $\kappa>0$ for all $c<j$, where
$\kappa$ is given by Eq. (LABEL:main-eq:kappa). This condition, however, still
allows for some freedom in the choice of $\alpha$, which can be used to reduce
the influence of the sub-leading orders to the Eq. (LABEL:main-eq:scaling).

These sub-leading orders are represented by the big-$\mathcal{O}$
contributions in Eqs. (S.7) and (S.8) and, especially for small supercells,
can change the behavior of certain bands from what would be expected based on
Eq. (LABEL:main-eq:scaling). Since we determine the confinement dimensionality
by analyzing the change of $\widetilde{V}^{\alpha}/\widetilde{E}$ with respect
to the reference supercell, the bands most susceptible to the sub-leading
order effects will be the bands where the smallest leading-order change is
expected. According to Eq. (LABEL:main-eq:scaling), these are the bands, for
which $\kappa$ is close to zero, i.e., for every $j$, the bands with $c=j$ and
$c=j-1$. Because we generally do not know the contribution of the sub-leading
orders to the overall scaling result, the most reasonable choice is to
maximize the leading-order movement, i.e., that of the $c=j$ bands downwards
and that of the $c=j-1$ bands upwards, at the same time.

Since $\kappa\colonequals-(\alpha+1)c+D$ is a linear function of $c$, the
desired leading-order movements will be maximal if $\kappa=0$ at the middle
point, i.e., for $c=j-1/2$. Thus, for each $0<j\leq D$, we solve the equation

$-(\alpha+1)(j-\frac{1}{2})+D=0.$ (S.10)

The solution to (S.10) is

$\alpha=\frac{2D}{2j-1}-1.$ (S.11)

By choosing $\alpha$ according to (S.11) for each $D,j$, we ensure the
effective minimization of the sub-leading order contributions to our scaling
analysis and thus the best precision of our method. The powers $\alpha$
obtained this way for each $j$ for $D=1,2,3$ are tabulated in Table
LABEL:main-table:kappa.

##  VI Electronic band computation

For our analysis of the 2D hexagonal BN layer with N vacancy, we calculated
the band structure and charge densities using the density functional theory
[10] implemented in the VASP code v6.1.1 [11]. We use the PBE exchange-
correlation functional [12] and represent the electron-ion interaction via the
PAW potentials [13], specifically via the B and N potential files recommended
by VASP. The cutoff energy was set to 550 eV.

We modeled the material using a 3D supercell with periodic boundary
conditions, starting each geometry optimization with the vacuum layer of 12.56
$\mathrm{\mathring{A}}$ perpendicular to the material layer to eliminate the
effects of periodic boundary conditions in this direction. For geometry
optimization and computation of the charge densities, we used $11\times 11$
$\Gamma$-centered $k$-point grid. Geometry optimization was performed via the
conjugate gradient algorithm until the residual atomic forces were smaller
than 10 $\mathrm{meV/\mathring{A}}$. For the computation of the band structure
in Fig. LABEL:main-fig:DFTconfinement(a), we interpolated each segment of the
high-symmetry path with 13 $k$-points, in addition to the two high-symmetry
points.

For completeness, we include the plot for $j=1$, which was omitted in the main
text, in Fig. S-2.

Figure S-2: Confinement analysis applied to electronic waves in a 2D hexagonal
BN with nitrogen deficiency, $j=1$. Zero energy is set at the Fermi level. (a)
Band structure of the $N=5$ supercell, flipped to have an energy abscissa.
Different bands are shown with different color for clarity. (b) Scaling
analysis of confinement. Every band is represented by a point. These results
confirm the three point-confined bands around the Fermi level. However,
compared to Fig. LABEL:main-fig:DFTconfinement, several other bands appear
slightly under the reference line. We attribute this effect to be most likely
caused by the sub-leading order contributions in the scaling and would be
remedied for larger supercell sizes. This conclusion is enhanced by the fact
that the standard deviation for these bands (green area) is higher than their
separation from the reference line.

These results confirm the three point-confined bands ($c=2$ in a $D=2$ system)
around the Fermi level. At the same time, some other bands of the $N=5$
supercell seem to be slightly under the reference values. Naively, this could
indicate $c=1$ confinement dimensionality for those bands, however, in this
case we attribute this effect to be most likely caused by the sub-leading
order effects in the scaling. This would be remedied by using larger supercell
sizes $N$ and $N_{0}$. This conclusion is confirmed by looking at the standard
deviation of the $c\neq 2$ bands for the $N=5$ supercell (green area in Fig.
S-2(b)), which is higher than the separation of the questionable bands from
their reference line, in each case. Fig. S-2 exemplifies that our framework
can offer rigorous and cost-effective scaling analysis, but, especially for
small supercell sizes, its results must still be critically evaluated to avoid
the sub-leading order effects of scaling.

##  VII Photonic band computation

In the main text we analyze light confinement in a 3D inverse woodpile
photonic band gap superlattice with two proximate linear defects. The photonic
band gap crystal consists of bulk silicon ($\varepsilon=12.1$, see Ref. [14]),
in which nanopores of radius $R$ filled with air ($\varepsilon=1$) are etched
[15, 16]. We model the unperturbed crystal using a tetragonal unit cell with
lattice parameters $a$ in the $y$ direction and $b=\nicefrac{{a}}{{\sqrt{2}}}$
in the $x$ and $z$ directions. This unperturbed structure is depicted in Fig.
S-3(a).

Figure S-3: (a) Structure of an unperturbed inverse woodpile photonic crystal.
We use a tetragonal unit cell with lattice parameters $a$ in the $y$ direction
and $b=a/\sqrt{2}$ in the $x$ and $z$ directions. The pore radius is denoted
by $R$. The figure shows an $N=2$ supercell, with one unit cell designated by
the dashed border. (b) Design of the defect. The radius $R^{\prime}$ of two
proximate defect pores (shown in green) is altered. At their intersection a
region with excess of one material is created that serves as a point defect
(red glow).

A defect can be incorporated in the structure by altering the radius
$R^{\prime}$ of two proximate nanopores, as illustrated in Fig. S-3(b). In our
example from the main text we use the unperturbed pores of the radius
$R=0.24a$ and the defect pores with the increased radius $R^{\prime}=0.35a>R$.
This results in creation of a region with excess air at the intersection of
the two defect pores, serving as a point defect. We include one defect of this
type per supercell. Such a defect results in splitting of the defect bands
from the bottom edge of the band gap, in analogy with acceptor-doped lattices
in solid state physics.

We investigate the confinement of $N=4$ supercell and use $N_{0}=2$ as a
reference. The spectrum of the investigated cavity superlattice is computed by
the plane-wave expansion method using the well-known MPB code [17] and is
depicted in the main text in Fig. LABEL:main-fig:photonicconfinement(a). For
the band structure, we follow the high-symmetry path of the tetragonal unit
cell according to Ref. [18]. In order to use unified notation, we re-labeled
our high-symmetry points and axes in Fig. LABEL:main-
fig:photonicconfinement(a) according to Ref. [18], which assumes that the two
equal lattice constants are in the $x$ and $y$ directions.

We perform scaling analysis to identify the confinement dimensionality $c$ of
each band, as described by our theory introduced in this Letter. Each step of
the analysis for this $D=3$-dimensional system is described in Fig. S-4. We
choose the auxiliary power $\alpha$ based on Table LABEL:main-table:kappa.

Figure S-4: Step-by-step confinement analysis of bands for a 3D inverse
woodpile photonic crystal with two proximate line defects. Every band is
represented by a point. The values of the auxiliary power $\alpha$ for each
$j$ are (a) $j=3,\alpha=1/5$, (b) $j=2,\alpha=1$, (c) $j=1,\alpha=5$, as per
Table LABEL:main-table:kappa.

Fig. S-4(a) describes $j=3$, with $\alpha=\nicefrac{{1}}{{5}}$. Out of the
whole spectrum of unidentified bands, three bands around $\tilde{\omega}=0.6$
decrease in value of $\log(V^{\nicefrac{{1}}{{5}}}/E)$ with respect to the
reference supercell. We thus identify these bands as having $c=3$ confinement
dimensionality. This is in stark contrast with the statement made in Ref. [19]
that no point-confined bands exist in acceptor-like 3D inverse woodpile
crystals.

Fig. S-4(b) describes $j=2$, with $\alpha=1$. In this case, several additional
bands between roughly $\tilde{\omega}=0.5$ and $\tilde{\omega}=0.6$ drop below
the values of the reference supercell. These newly separated bands are
therefore identified as having $c=2$ confinement dimensionality. This marks
the first time such $c=2$ bands are distinguished from fully extended bands in
a 3D inverse woodpile photonic structure.

Fig. S-4(c) describes $j=1$, with $\alpha=5$. There are several more bands
dropping below the reference line. According to our framework, this would mean
that they are plane-confined, i.e., having $c=1$. Nevertheless, we know that
our defect geometry is linear and therefore cannot support plane-confined
bands. These bands must thus either have $c=2$ or $c=0$.

The misidentification for bands that should have $c=2$ is due to the fact that
the size of the (reference) supercell is still below the localization length
in the plane where the linear defects are positioned. However, the scaling Eq.
(LABEL:main-eq:scaling) holds only for supercell sizes larger than the
localization length. Therefore it is understandable that our method "thinks"
that these bands are plane-confined. Simply put, for the given supercell size,
the information about confinement of these bands is not yet available in the
calculated energy density and thus cannot be resolved by any means, unless
implementing some additional knowledge, for example about the defect geometry.

The misidentified bands that should have $c=0$ are caused by the fact that for
the given supercell size the effects of the sub-leading order terms on Eq.
(LABEL:main-eq:scaling) are still too large. This is an issue of the
convergence of our technique and it would also be corrected for larger
supercell sizes. In this case, it could be possible to devise some alternative
way to improve the convergence since, unlike in the case of the previous
paragraph, the information about the confinement may be available in the
calculated energy density data - this remains, however, beyond the scope of
this Letter.

Lastly, all the bands remaining above the reference values in Fig. S-4(c) are
identified as extended, i.e., with $c=0$.

Our method is especially powerful for experimentally interesting supercells of
moderate size. These supercell sizes are, however, in a range where possible
inaccuracies due to finite-size scaling can occur and therefore, naturally,
the results need to be critically evaluated. The presented photonic example
further illustrates the strength and novelty of our method. Our framework has
enabled the discovery of 3D-confined acceptor-like bands previously thought
non-existent and, for the first time ever in an inverse woodpile structure,
distinguished 2D-confined bands from the extended ones. It is therefore clear
that our technique enables direct access to previously inaccessible
confinement information.

## References

  * Ku  _et al._ [2010] W. Ku, T. Berlijn, and C.-C. Lee, Phys. Rev. Lett. 104, 216401 (2010).
  * Mayo  _et al._ [2020] S. G. Mayo, F. Yndurain, and J. M. Soler, J. Phys. Condens. Matter 32, 205902 (2020).
  * Aoki [1983] H. Aoki, J. Phys. C: Solid State Phys. 16, L205 (1983).
  * Janssen [1998] M. Janssen, Phys. Rep. 295, 1 (1998).
  * Faez  _et al._ [2009] S. Faez, A. Strybulevych, J. H. Page, A. Lagendijk, and B. A. van Tiggelen, Phys. Rev. Lett. 103, 155703 (2009).
  * Hurley and Rickard [2009] N. Hurley and S. Rickard, IEEE Trans. Inform. Theory 55, 4723 (2009).
  * El-Dardiry  _et al._ [2012] R. G. S. El-Dardiry, S. Faez, and A. Lagendijk, Phys. Rev. B 86, 125132 (2012).
  * Arceri and Corwin [2020] F. Arceri and E. I. Corwin, Phys. Rev. Lett. 124, 238002 (2020).
  * Pashartis and Rubel [2017] C. Pashartis and O. Rubel, Phys. Rev. B 96, 155209 (2017).
  * Kratzer and Neugebauer [2019] P. Kratzer and J. Neugebauer, Front. Chem. 7, 106 (2019).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Computational Materials Science 6, 15 (1996).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Hillebrand  _et al._ [2003] R. Hillebrand, S. Senz, W. Hergert, and U. Gösele, J. Appl. Phys. 94, 2758 (2003).
  * Ho  _et al._ [1994] K. M. Ho, C. T. Chan, C. M. Soukoulis, R. Biswas, and M. Sigalas, Solid State Commun. 89, 413 (1994).
  * Leistikow  _et al._ [2011] M. D. Leistikow, A. P. Mosk, E. Yeganegi, S. R. Huisman, A. Lagendijk, and W. L. Vos, Phys. Rev. Lett. 107, 193903 (2011).
  * Johnson and Joannopoulos [2001] S. Johnson and J. Joannopoulos, Opt. Express 8, 173 (2001).
  * Setyawan and Curtarolo [2010] W. Setyawan and S. Curtarolo, Comput. Mater. Sci. 49, 299 (2010).
  * Woldering  _et al._ [2014] L. A. Woldering, A. P. Mosk, and W. L. Vos, Phys. Rev. B 90, 115140 (2014).