Scaling theory of wave confinement in classical and quantum periodic systems

Functional defects in periodic media confine waves - acoustic, electromagnetic, electronic, spin, etc. - in various dimensions, depending on the structure of the defect. While defects are usually modelled by a superlattice with a typical band-structure representation of energy levels, determining the confinement associated with a given band is highly non-trivial and no analytical method is known to date. Therefore, we propose a rigorous method to classify the dimensionality of the confinement. Starting from the confinement energy and the mode volume, we use finite-size scaling to find that ratios of these quantities to certain powers yield the confinement dimensionality of each band. This classification has negligible additional computational costs compared to a band structure calculation and is valid for any type of wave in both quantum and classical regimes, and any dimension. In the quantum case, we illustrate our method on electronic confinement in 2D hexagonal BN with a nitrogen vacancy, which confirms the previous results. In the classical case, we study a three-dimensional photonic band gap cavity superlattice, where we identify novel acceptor-like behavior.

Functional defects in periodic media confine waves -acoustic, electromagnetic, electronic, spin, etc. -in various dimensions, depending on the structure of the defect. While defects are usually modelled by a superlattice with a typical band-structure representation of energy levels, determining the confinement associated with a given band is highly non-trivial and no analytical method is known to date. Therefore, we propose a rigorous method to classify the dimensionality of the confinement. Starting from the confinement energy and the mode volume, we use finite-size scaling to find that ratios of these quantities to certain powers yield the confinement dimensionality of each band. This classification has negligible additional computational costs compared to a band structure calculation and is valid for any type of wave in both quantum and classical regimes, and any dimension. In the quantum case, we illustrate our method on electronic confinement in 2D hexagonal BN with a nitrogen vacancy, which confirms the previous results. In the classical case, we study a threedimensional photonic band gap cavity superlattice, where we identify novel acceptor-like behavior.
The analysis of spatial concentration of energy in physical systems is in photonics traditionally done via the mode volume [18,[36][37][38][39] or equivalently, in condensed matter physics, via the participation ratio [16,23,40]. Bands with small mode volume or participation ratio are considered to be confined while, conversely, bands with large mode volume or participation ratio are taken to be extended [16,18,19,[22][23][24]41]. However, the notion of what specifically is 'large' and 'small' is in each case determined rather subjectively and there are no rigorous boundaries imposed by nature.
An alternative method to analyze the wave confinement, called multifractality analysis [42][43][44], is based on the participation ratio scaling in the limit of infinitely large supercells. Unfortunately, this approach requires impractically large supercells to yield reliable results, further compounded by its inability to deal with band folding (see Supplemental material).
In this Letter, we present a rigorous method to determine confinement of waves in periodic structures with defects, based on finite-size scaling [45][46][47]. Rather than extending the supercells towards infinity, our method determines confinement by looking at supercells of finite size. As a consequence, our approach requires only minimal computational overhead and its results are directly relevant for experimentally interesting finitely-large systems. This technique can be viewed as a more accessible and practical extension of the multifractality concept, also suitable for automated classification.
Superimposing a defect lattice on an unperturbed crystal lattice gives rise to a defect superlattice [6,22,23,[48][49][50]. A unit cell of such a superlattice, containing N D unit cells, where D is the system dimension, will be referred to as a supercell of linear size N . For simplicity, we keep N uniform for all directions.
Based on the geometric dimensionality d of the defects, wave bands with different confinement dimensionality c arise in the spectrum, as depicted in Fig. 1 for an N = 5 supercell in D = 3. The confinement dimensionality c of a wave determines the number of dimensions in which the wave is confined and mathematically corresponds to the codimension of the defect dimensionality: c = D − d.
We note that introducing defects in every unit cell of the supercell is equivalent to a new periodic structure with no defects at all (d = 3) and supports only extended, i.e., unconfined (c = 0), wave bands. Fig. 1 represents highly idealized examples of simple defects. Real structures often contain multiple defects of various geometries oriented along different directions. Such complex defect configurations allow for several different confinement dimensionalities c and the problem of assigning the correct c to a given band thus becomes highly nontrivial. In order to understand the confinement in these real structures, a systematic method of identification and classification of confined bands is needed. (a) Four supercells, each containing a point defect akin to a cavity [11], trapping waves in all three dimensions, d = 0, c = 3. (b) Supercell with a linear defect, analogous to an optical waveguide or fiber, d = 1, c = 2. (c) Supercell with a planar defect, analogous to a 2D electron gas [51], d = 2, c = 1.
To analyze the confinement of waves, one must first determine the physical quantity corresponding to the notion of wave confinement for the given type of physical wave. Without loss of generality, one can make this a real, non-negative quantity. For example, in electromagnetic systems this quantity could be the energy density, in electronic systems it could be the charge density, and similarly for other waves. To keep the description of our method general, in the following we will assume that this quantity has been identified and that it can be calculated for every band of interest. We will denote this spatially dependent confinement quantity by W (x) and it will serve as the starting point of our analysis.
For scientifically interesting structures of moderate supercell sizes, the number of bands to be analyzed often reaches several hundreds (see, e.g., [18,19]). It is therefore clearly impractical, borderline unfeasible, to analyze the confinement by visually inspecting the spatial distribution of W (x) for each band. Such an approach is also prone to human error and, due to its qualitative nature, hard to automate. There is thus a clear need for a more quantitative approach towards the wave confinement analysis.
We employ W (x) to define the two key quantities for our scaling analysis. First is the mode volume where the integration is over the supercell volume V S . We note that in literature the definitions of the mode volume may differ for various applications [18,[36][37][38][39]. This difference is usually manifested in modification of the integration volume in Eq. (1). In our case of a superlattice, it seems physically obvious to choose the supercell volume V S as the integration domain. Moreover, as we show below, irrespective of the physical meaning of the quantity V M , integrating over V S is crucial for obtaining the correct scaling behavior in our theory. As an alternative quantity to the mode volume, the participation ratio [16,40] can be used, which also involves integration over V S . For details on this alternative, see Supplemental material.
The mode volume intuitively corresponds to the volume in which the wave is confined. On the other hand, it is important to also consider the problem from the converse point of view: How much energy of the wave is stored in a certain volume? In order to quantify this notion, we define the confinement energy with the integration volume V C taken as the volume of one unit cell, centered around the defect. An analogous quantity [52] has been introduced by Ref. [53] to analyze confinement of surface acoustic waves. For practical purposes, we employ normalized quantities, in the following: Here, E S denotes the total amount of energy in the supercell. It is clear that for a wave band confined within a cavity we expect simultaneously low V and high E, while the opposite (high V , low E) is expected for a fully extended band. Nevertheless, there are no natural thresholds on how low V and how high E should be for a band to be clearly identified as confined. To overcome this ambiguity, we aim to analyze the behaviour of V and E with respect to the variation of the supercell size N , instead of their values themselves. This technique is known as finite-size scaling [45][46][47].
Our scaling argument is illustrated with the following didactic 1D model in Fig. 2. Let us consider a supercell with N = 4 with a cavity in its center. For a specific band with certain V and E, we are interested in how adding more unit cells to the supercell boundary changes these values, and how these changes differ for a confined band as opposed to an extended band. Therefore, we increase the supercell size to N = 6. While scaling, we keep the total amount of energy within the supercell constant, i.e., E S = const.
A confined band is depicted in Fig. 2(a) with its energy density decaying away from the cavity in the supercell center. The size of the cavity to which the wave is confined and thus also the mode volume V M remain constant while scaling. However, the supercell volume V S has increased proportionally to N . Hence, the normalized quantity V decreases for the confined band as N −1 . Moreover, since the wave is confined within the cavity, it only 'feels' that additional unit cells have been added through its decaying tail. This response clearly approaches zero as N → ∞ and, thus, in this limit E = const.
An extended band is illustrated by a constant energy density distribution throughout the whole supercell in Fig. 2(b). In this case, the behavior described above is the converse: The volume occupied by the wave in the supercell grows proportionally to the supercell volume, resulting in the scaling V = const. The energy density extends homogeneously throughout the whole supercell and will further spread into the added unit cells as N is increased. This leakage decreases the amount of confine- It is straightforward to extend the above scaling analysis to a D-dimensional system. A band with a given confinement dimensionality 0 ≤ c ≤ D obeys the following scaling relations, in the limit of N → ∞: Here, A and B are constants independent of N . For derivation of Eqs. (4), see Supplemental material. One can directly calculate V and E from Eqs. (1)-(3) for each band, but because the constants A, B are not known a priori, Eqs. (4) cannot be easily inverted to obtain c. To accurately obtain c from our scaling relations, we combine the two equations in (4) into a ratio of the normalized mode volume raised to a judiciously chosen power α > 0 and the normalized confinement energy: where the critical exponent κ is given by The equation (5) strictly holds only in the limit of very large N , with sub-leading order terms in N present for finite supercell sizes. However, in Eqs. (5) and (6), we have introduced the auxiliary power α, which can be utilized to effectively minimize the contribution of these sub-leading terms and thus identify the confinement of bands already for very small supercells. Specifically, by suitably choosing the power α, we can adjust the value of κ so that it is negative for bands with certain c and positive for the other bands. The value of κ, in turn, influences the scaling behavior of these bands as per Eq. (5).
Our technique is illustrated in the flow diagram in Fig. 3. Investigating the confinement in a supercell of size N , for every 0 < j ≤ D, we have chosen α so that κ < 0 for all c ≥ j and κ > 0 for all c < j. We are now able to distinguish between the bands with c < j and c ≥ j by simply tracking the behavior of V α / E as the supercell grows from a smaller reference size N 0 to the investigated size N . The bands with negative κ will move downwards in the graph, while the bands with positive κ will shift upwards, in accordance with the Eq. (5). By varying j over all the integer values 0 < j ≤ D, this approach allows us to completely classify all wave bands in the spectrum based on their confinement dimensionality.
The unique values of α to maximize the accuracy of our technique for D = 1, 2, 3 are listed in Table I. For details on their calculation, see the Supplemental material. In case of the 1D example from Fig. 2, if we choose α = 1, corresponding to D = 1 and j = 1 from Table I, upon plottingṼ /Ẽ versus the supercell size N we will observe that the c = 1 confined band moves down, while the c = 0 extended band moves up in the flow diagram, analogously to Fig. 3.
We now demonstrate our method on a 2D hexagonal BN with a nitrogen vacancy representing a point-like defect [54]. Specifically, we investigate electronic confinement in a supercell of size N = 5. The band structure and charge densities of this quantum system were calculated using density functional theory [55] implemented in the VASP code v6.1.1 [56]. For details see Supplemental material.
The band structure of the system is shown in Fig. 4(a). Since the defect in our D = 2 system has a point geometry, we only expect two different confinement dimensionalities to appear: the point-confined (c = 2) and the extended (c = 0) bands. To distinguish between these, we choose α = 1 /3 from Table I  Flow diagram to identify bands confined in c ≥ j dimensions in a D-dimensional medium: The ratio V α / E versus the system size N . The auxiliary power α is chosen so that κ < 0 for all c ≥ j and κ > 0 for all c < j. The sign of κ determines if V α / E increases or decreases with growing supercell size. Fig. 4(b) clearly shows that the majority of bands move upwards, identifying them as extended (c = 0) bands according to our framework. Additionally, three bands near zero energy move downwards and thus correspond to point-confined (c = 2) bands.
In our second example, we demonstrate our technique on light confinement in an N = 4 supercell of a 3D inverse woodpile photonic crystal with two proximate line defects [18,19,57,58]. The crossing point of these line defects represents a point defect. We specifically study the acceptor-like structure investigated by Ref. [18]. It is clear that the defects in this case have more complicated geometry than in the previous example: The system sus- bands and extended (c = 0) bands. The band structure and the energy densities were calculated with the MPB code [59]. For details on the structure and analysis, see Supplemental material. Fig. 5 depicts the results of our analysis. We unambiguously identify the confinement dimensionality c of most bands. Additionally, several bands, mostly near ω ≈ 0.5, appear to be plane-confined (c = 1). We know, however, that our defect geometry is linear and thus it does not support plane-confined bands. We attribute this discrepancy to the small size of the studied supercell. More specifically, these bands either have c = 0, with the sub-leading order contributions to Eq. (5) taking over for the studied supercell sizes, or they have c = 2, with localization length in one of the dimensions still being larger than the reference supercell size N 0 , thus appearing as only being confined in one dimension. Nevertheless, our analysis is the first that distinguishes c = 2 bands from those with c = 0 in a 3D photonic superlattice.
We also find three bands (red) with c = 3. Visual inspection confirms the point-confined character of the two bands atω ≈ 0.62. The inspection also shows that the band atω ≈ 0.56 is confined in at least two dimensions, however it is visually not clear if the confinement is c = 2 or c = 3. Nevertheless, our analysis clearly disproves the statement in Ref. [18] that no point-confined bands exist in acceptor-like inverse woodpile crystals.
In this Letter we described a systematic scaling theory to analyze wave confinement in periodic media with functional defects, applicable to any type of physical waves. Already in one of the studied examples our technique uncovers physically new results, thus showcasing its power. Our method is directly applicable to actively researched periodic structures and for optimization algorithms aiming to minimize or maximize specific types of wave confinement.
We thank Geert Brocks and Menno Bokdam for access to the CMS cluster with VASP code to perform DFT calculations. We also thank Geert Brocks for recommending hexagonal BN as a suitable material for illustrating our framework. Furthermore, we thank Sjoerd Hack for the introduction into the MPB software, and Lars Corbijn van Willenswaard and Manashee Adhikary for stimulating discussions on the physics of waves. This re-search is supported by the Shell-NWO/FOM programme "Computational Sciences for Energy Research" (CSER) and the MESA+ Institute for Nanotechnology, Applied Nanophotonics (ANP).

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 V α / 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. (4). 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 α, 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 C in (2) we choose the volume V C to be centered around the cavity, where the confinement is expected. Nevertheless, this is not a strict requirement. Choosing the volume V 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 C is selected, to keep this volume constant throughout the scaling, so that the exponent c − D appears in the scaling relation for E in (4), since this is crucial for being able to tune the values of the critical exponent κ 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. (1) and (2) must be performed over its whole three-dimensional volume, i.e., the volumes V S and V C will be three-dimensional. Nevertheless, since the structure of the material is still twodimensional, 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 ≤ 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 α in case of fractional d or arbitrarily high D. These possibilities further emphasize the generality of our scaling theory.

SPARSITY MEASURES: ANALOGIES OF MODE VOLUME AND PARTICIPATION RATIO
In our confinement analysis, we employ the mode volume V M as one of two key quantities, defined in Eq. (1). This is, however, not the only possible choice, and a whole class of quantities can be substituted for V 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 , . . . , m M } ∈ R M , is a pq-mean: Ref. [6] has shown that, for p < q < ∞, µ p,q is indeed a viable measure of sparsity. There is clear correspondence between the inequality of W (x) over the volume V S for a given band and its confinement: For a confined band, W (x) is very high in certain parts of V S but low everywhere else, while for an extended band W (x) is approximately equally distributed over the whole V S . To see this quantitatively, we write arXiv:2205.00514v1 [cond-mat.dis-nn] 1 May 2022 the straightforward continuous generalization of (S.1) to quantify the inequality of W (x) over the volume V S as Our normalized mode volume V , as defined by Eq. (3), turns out to be a continuous version of the pq-mean with p = 1 and q = ∞, i.e., V = U 1,∞ . The fact that q = ∞ means that 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 V and our method will still work properly, albeit possibly with different scaling power than we obtained for V in Eq. (4). Specifically, substituting p = 1, q = 2 yields The quantity U 1,2 corresponds to the square root of normalized participation ratio: where the participation ratio P [7-9] is given by 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 2 1,2 instead of U 1,2 . This illustrates that, in our approach, instead of normalized mode volume V , normalized participation ratio can be used as well, similarly to any other pq-mean U p,q .

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 For some bands the plotted values change significantly when scaling from N0 = 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 V α / E. use the N 0 = 2 supercell as a reference. We scale either the mode volume V , in Fig. S-1(a), or the confinement energy E, in Fig. S-1(b).
From Fig. S-1(a), we observe that for some bands the value of V decreases significantly as the supercell grows and for some it stays roughly constant. One would immediately expect the bands with constant V to be extended, according to Eq. (4). 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ω ≈ 0.5.
Similarly, in Fig. S-1(b), for some bands the value of 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 V and 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 V and E.

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 ≤ 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 (1), for such a wave is given by where A is a constant independent of N . In fact, A is actually a functional depending on the permittivity distribution (x), the wave-vector k and the frequency ω.
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 S = N D , we obtain Similarly, the normalized confinement energy, defined by (2) and (4), behaves as where B is a constant independent of N . Analogously to A, B is in fact a functional depending on the permittivity distribution (x), the wave-vector k and the frequency ω. 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 O(N c−D−1 ). We combine the Eqs. (S.7) and (S.8) and, upon neglecting the big-O contributions, we obtain for a wave with the confinement dimensionality c in the limit N → ∞: where C = A α /B. Eq. (S.9) represents exactly our scaling relation in Eq. (5).

DETERMINING THE AUXILIARY POWERS α
In our method, we are able to distinguish between the bands with c < j and c ≥ j by tracking the behavior of V α / E according to Eq. (5), 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 ≤ 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 α has to be chosen so that κ < 0 for all c ≥ j and κ > 0 for all c < j, where κ is given by Eq. (6). This condition, however, still allows for some freedom in the choice of α, which can be used to reduce the influence of the sub-leading orders to the Eq. (5).
These sub-leading orders are represented by the big-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. (5). Since we determine the confinement dimensionality by analyzing the change of V α / 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. (5), these are the bands, for which κ 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 κ := −(α + 1)c + D is a linear function of c, the desired leading-order movements will be maximal if κ = 0 at the middle point, i.e., for c = j − 1/2. Thus, for each 0 < j ≤ D, we solve the equation The solution to (S.10) is By choosing α 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 α obtained this way for each j for D = 1, 2, 3 are tabulated in Table I.

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  Fig. 4, 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.
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Å 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 × 11 Γ-centered k-point grid. Geometry optimization was performed via the conjugate gradient algorithm until the residual atomic forces were smaller than 10 meV/Å. For the computation of the band structure in Fig. 4(a), we interpolated each segment of the highsymmetry 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. 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 = 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 subleading order effects of scaling.

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 (ε = 12.1, see Ref. [14]), in which nanopores of radius R filled with air (ε = 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 = a / √ 2 in the x and z directions. This unperturbed structure is depicted in Fig. S-3(a). A defect can be incorporated in the structure by altering the radius R 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 = 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. 5(a). For the 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 α for each j are (a) j = 3, α = 1/5, (b) j = 2, α = 1, (c) j = 1, α = 5, as per Table I. 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. 5(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 α based on Table I. Fig. S-4(a) describes j = 3, with α = 1 /5. Out of the whole spectrum of unidentified bands, three bands aroundω = 0.6 decrease in value of log V 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 α = 1. In this case, several additional bands between roughlyω = 0.5 and ω = 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 α = 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. (5) 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. (5) 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.