Transition from light diffusion to localization in three-dimensional amorphous dielectric networks near the band edge

Localization of light is the photon analog of electron localization in disordered lattices, for whose discovery Anderson received the Nobel prize in 1977. The question about its existence in open three-dimensional materials has eluded an experimental and full theoretical verification for decades. Here we study numerically electromagnetic vector wave transmittance through realistic digital representations of hyperuniform dielectric networks, a new class of highly correlated but disordered photonic band gap materials. We identify the evanescent decay of the transmitted power in the gap and diffusive transport far from the gap. Near the gap, we find that transport sets off diffusive but, with increasing slab thickness, crosses over gradually to a faster decay, signaling localization. We show that we can describe the transition to localization at the mobility edge using the self-consistent theory of localization based on the concept of a position-dependent diffusion coefficient.

B andgap formation and strong Anderson localization (SAL) of classical waves are both considered general wave phenomena where the mechanism leading to the exponential attenuation of wave transport can be understood in terms of interference of scattered waves. In a periodically repeating environment, scattering of waves from Bragg planes can be associated with the opening of a photonic bandgap (PBG) [1][2][3] .
The SAL mechanism is usually explained by the constructive interference of multiply scattered waves propagating along timereversed loops, which increases the return probability and eventually leads to a breakdown of wave diffusion [4][5][6] . In contrast to PBG formation, the transition to SAL in disordered media strongly depends on dimensionality. In one and two dimensions, there are no truly extended states, and waves in disordered media are always localized for sufficiently large systems sizes [7][8][9][10][11] . Only in three dimensions, wave localization shows a phase transition, and localized states appear below the 'mobility edge'. The threshold for localization can be estimated by the Ioffe-Regel criterion k·ℓ~1, where ℓ denotes the transport mean free path and k the wavenumber [12][13][14] . The Anderson transition in three dimensions is a fascinating phenomenon that until now has eluded a full theoretical and experimental understanding. It's relevance is not restricted to transport of photons or electrons, but it can also be applied in acoustics, or any kind of coherent wave propagation 15 . Moreover, a full understanding and control over the wave fields in complex media offer a plethora of opportunities for applications in imaging, sensing, and photonics 6 . However, an experimental observation of SAL for electromagnetic vector waves in threedimensional systems has not yet been achieved, even though several claims to its existence were made [16][17][18] but soon after were put in question 19,20 and later refuted 21 . At about the same time it was found theoretically that SAL of electromagnetic vector waves is absent in random ensembles of point scatterers, irrespective of their scattering strength 22 . Thus, in a 2016 perspective article, Skipetrov and Page 21 declare a "Red light for Anderson localization" of light and Maret et al. 23 ask "Can 3D light localization be reached in 'white paint'" at all?
The advent of amorphous PBG materials over the last decade [24][25][26][27][28][29] has opened a new pathway toward the design of strongly photonic dielectric materials. It has also prompted fundamental questions concerning the relationship between the Anderson localization transition and the transition to a full bandgap. Recently, we proposed a transport phase diagram for twodimensional hyperuniform disordered PBG materials with a SAL regime near the PBG 30 and conjectured that these findings could be generalized to three dimensions [24][25][26][27]30,31 . Early pioneering work by John suggested the possibility of finding SAL in disordered crystalline structures near the bandgap 2,32,33 , an idea supported more recently by the numerical studies of Conti and Fratalocci 33 . Imagawa et al. reported an increase of the inverse participation ratio near the bandgap of an amorphous diamond structure, which might indicate the presence of localized states 34 . This previous work provides a rationale for the existence of a SAL regime in the vicinity of a bandgap of an amorphous photonic material, which we are going to investigate in our work.
To this end, we study the transport properties of electromagnetic vector waves in realistic digital representations of threedimensional hyperuniform silicon networks numerically. We find evidence for the anomalous light transport near the band edge, signaling the onset of Anderson localization at a mobility edge and a broad frequency window where light is localized before the bandgap fully develops.

Results
Strong Anderson localization of waves. The discussion of what defines SAL is very detailed and rich, and for a comprehensive review, we refer to the literature 5,13,21,32,[35][36][37][38] . A working definition for finite-sized systems, proposed by Cherroret and Skipetrov, is that 'SAL is an interference wave phenomenon in a medium of a finite size that would give rise to truly localized states if the medium were extended to infinity' 38 . The transition to SAL in three dimensions is usually described in the framework of the self-consistent (SC) theory 39 . It treats localization by introducing a position r ! dependent wave diffusion coefficient Dð r ! Þ, which decays to zero deep inside the medium 5,38 , due to an increased return probability of multiple scattering paths. At the mobility edge, for a semi-infinite medium, one finds the simple algebraic forms DðzÞ ¼ Dð0Þ 1þz=ξ c and in the localized regime DðzÞ ¼ Dð0Þ expðÀ2z=ξÞ where z denotes the distance from the surface. The localization length ξ becomes finite at the critical point while D(z) < D B already when approaching the mobility edge, where D B = cℓ/3 denotes the standard Boltzmann diffusion constant. D(0) continues to drop gradually as the localization threshold is crossed 40 . For slabs of finite thickness L, the transmittance T(L) shows two distinct regimes: initially, it decays diffusive as~L −1 which is followed by an exponential decaỹ e −L/ξ . Precisely at the transition, SC theory predicts a critical power-law scaling T~L −2 instead of the exponential decay.
We point out that the (SC) theory is an approximate theory, and different variants have been proposed in the literature 5,38,41 . As discussed in ref. 38 , in its original form, it is strictly valid only for k ⋅ ℓ ≫ 1. The fact it can describe certain phenomena at the mobility edge and in the localized regime is somewhat fortuitous and not fully understood. However, progress toward a better microscopic understanding has been reported recently 38,42 . Moreover, almost all theoretical and numerical studies were carried out assuming scalar waves and point scatterers randomly distributed in space, i.e., using white-noise Gaussian statistics 14 . It is not self-evident that the conventional SC theory 40 , can explain the transport of vector waves in spatially correlated, densely filled dielectric materials with a bandgap. First of all, the transport ℓ and the scattering mean free path ℓ s are not the same for scatterers of size~λ, and the wave is propagating in some effective medium with a wavenumber k eff 11,30 . If and how this affects the predictions by SC theory is currently not known. As a consequence of the approximate nature of the SC theory, in particular, when applied to realistic representations of dielectric materials, we must assume that there is no perfect one-to-one relationship between the macroscopic transport properties, such as the localization length ξ and microscopic quantities like k = 2π/λ and ℓ. We, therefore, denote with (kl) the localization parameter, which we assume is similar and proportional but not necessarily identical to k⋅ℓ. For simplicity, we also assume (kl) c ≡ 1 14 , and we have tested that using slightly different values for the localization threshold does not significantly affect our findings.
In summary, in our study the localization parameter (kl) sets the macroscopic properties in the SC theory, such as ξ/ℓ = 6(kl) 2 / (1 − (kl) 4 ) for (kl) < 1 and D(0), with D(0) = D B (1 − (kl) 2 ) for (kl) ≫ 1 40 . The relationship to the microscopic parameter ℓ is established via (kl) ∝ k⋅ℓ modulus a prefactor of order one that also takes account of the uncertainty with respect to ℓ/ℓ s and the effective wavenumber k eff ≳ k. We note that the predictions by SC theory are only meaningful in the limit ξ/ℓ ≫ 1. For example ξ/ℓ > 5 for ðklÞ 2 ½ 3 4 ; 1Þ.
Transport in a photonic bandgap. A photonic crystal with a full PBG displays a vanishing density of states (DOS) in the gap. Transport through a finite-sized slab is due to tunneling, characterized by an exponential attenuation with a decay length L B , called the 'Bragg length,' typically on the order of one unit cell, in high refractive index photonic crystals 43,44 . For frequencies outside the gap, or for specific directions in the presence of an incomplete gap, photonic crystals are transparent, before the onset of diffraction 43 . Perfect crystals show neither photon diffusion nor localization. In his celebrated 1987 paper, John suggested that three-dimensional photonic crystal lattices with moderate disorder may exhibit strong localization of photons 2 , an idea supported more recently by numerical studies 33 . Amorphous PBG materials are disordered but the spatial distribution of dielectric material is highly correlated, which can also lead to the opening of a full PBG. The genuine disorder, however, implies that these materials display strong scattering for frequencies outside the gap, or if the gap is incomplete for all frequencies, with the notable exception of transparency in stealthy hyperuniform materials in the long-wavelength limit 45,46 , which we do not address here. Strong scattering outside the gap opens the possibility for the existence of SAL transport regimes, even in the absence of defect states 30 .
Optical transport simulations and density of states (DOS). We study the transport of waves in three-dimensional hyperuniform silicon photonic network structures, refractive index n = 3.6, derived from the center positions of an assembly of 10,000 randomly close-packed spheres, diameter a, as described earlier 47,48 .
The design protocol consists of mapping the seed pattern into tetrahedrons by performing a Delaunay tessellation. Then, the centers of mass of the tetrahedrons are connected, resulting in a tetravalent network structure of interconnected rods with the desired structural properties. The diameter of the rods sets the space-filling fraction ϕ, Fig. 1a, c. The length a is the typical shortrange structural length scale of the network, which for a crystal would be the lattice constant. The seed point pattern is hyperuniform, but not stealthy, meaning that the isotropic structure factor vanishes asymptotically in the limit S(k) → 0 for k → 0 where k = 2π/λ denotes the wavenumber in vacuum. Practically identical network structures have been considered by Liew et al. in a study of the optical DOS 26 (see also Supplementary Figs. 2 and 3). They report a substantial depletion of the DOS, by more than two orders of magnitude, over a significant range of frequencies, indicating the presence of a bandgap for different values of ϕ. In the present study, we consider networks with a filling fraction of ϕ = 0.28, shown to display the most pronounced photonic properties in their study 26 .
For the optical transport simulations, we apply the finite differences time-domain (FDTD) approach, implemented by the MIT Electromagnetic Equation Propagation (MEEP) 49 . It is considered to be one of the most potent simulation techniques to study electromagnetic wave transport. In a single MEEP-simulation run, a broadband pulse of linearly polarized light, with an electric field vector parallel to one of the sides of the simulation box, is incident on the sample, as illustrated in Fig. 1a. We obtain the full spectrally resolved information about the optical transmittance T(a/λ, L), Fig. 1b. We present all spectra in terms of the reduced frequency ν 0 :¼ a=λ ¼ νa=c where λ, ν, and c denote the wavelength, frequency and speed of light  Fig. 1 Numerical simulations of optical transport properties. a Cross section of a three-dimensional FDTD (MEEP) simulation box containing a hyperuniform silicon (n = 3.6) network structure, thickness L = 6a. A light wave, linearly polarized along the x-axis, is launched on the leftside and propagates along the z-axis. The photonic network structure is terminated with perfectly matched layers (PML) at both sides of the box along the propagation axis. PMLs act as absorbers. The source (SRC) and detector (MON) are placed at a distance approximately~a from the sample, which is held in vacuum. Periodic boundary conditions are applied along x and y directions. b Triangles: transmittance spectrum T(a/λ) for a slab of thickness L = 18a for a filling fraction ϕ = 0.28. The optical transport data is compared to numerical calculations of the density of states (DOS) (squares). The bandgap-center frequency is ν 0 Gap ¼ a=λ Gap ¼ 0:478 and the width Δν is indicated by the shaded area. c Three-dimensional rendering of a hyperuniform network structure, edge length 6a and filling fraction ϕ = 0.28. The size of the structure used in the simulation is 18a × 18a × L with L ≤ 18a, which is repeated periodically in (x, y) direction to construct the slab geometry. d In the gap the transmittance decays exponentially and ln T collapses on a master curve when plotted in reduced units L/L B . L B ≤a denotes the Bragg length and it is found to be smallest near the gap-center frequency ν 0 Gap ¼ 0:478, see also  5 . We obtain results for slabs of thickness L = 0.3a-18a and average the simulations results over 6 (thick slabs) to 15 (thin slabs) independent configurations of the network structure. We note that the networks are structurally isotropic 47 and therefore, all incident polarization states on average lead to the same results.
To obtain the normalized photonic DOS we use the supercell method 3 , implemented in the open-source code MIT Photonic Bands (MPB) 50 , Fig. 1b. Due to computational limitations, we have generated equivalent, but smaller seed patterns with periodic boundary conditions applied. To this end, we are using a packing algorithm developed by Skoge et al. 51 .
Fitting procedure. The transmittance spectrum Tðν 0 ; LÞ depends on a number of a priori unknown parameters. Extracting all these parameters from a global fit to SC theory is not stable and prone to overfitting. To overcome this problem, we first determine parameters that do not scale with ν 0 , such that eventually, Tðν 0 ; LÞ predicted by SC theory, depends only on the one parameter (kl). Once this is achieved, the predictions by SC theory and diffusion theory can be compared without fitting bias.
We test the presence of three transport regimes: (1) evanescent transport ln ½Tðν 0 ; LÞ $ ÀL=L B in a PBG with a Bragg length L B ℓ.
(2) Diffusive transport with DðzÞ % const: and (3) SAL, as described by the SC theory with (kl) ≲ 1 and ξ ≫ ℓ. We proceed as follows. (1) We start by looking for the signature of an evanescent decay in the gap. (2) Next, we consider frequency intervals far below and above the gap, where the predictions for T(L) from classical diffusion theory are sufficient to describe the data. This is the case whenever L ≫ ℓ and DðzÞ % const: From this fit we extract the angle averaged reflection coefficient R. (3) Next we identify the position(s) of the mobility edge(s) ν 0 c , where kl ð Þ ¼ kl ð Þ c 1, by fitting the data with SC theory treating both (kl) and and ℓ as adjustable parameters. From this fit we find the anchor points ν 0 c where kl ð Þ ¼ 1. Together with k ¼ 2π=λ ¼ 2πν 0 =a, this provides us with the proportionality between k ⋅ ℓ and (kl). (4) Now, that we have fixed all other parameters, we will attempt to describe the entire data set for Tðν 0 ; LÞ by SC theory with only one adjustable parameter ðklÞ½ν 0 .
Density of states and tunneling through the gap. Our numerical calculations of the DOS reveal a full bandgap for frequencies in the interval ν 0 2 ½0:47; 0:49, Fig. 1b (see also Supplementary Methods). Our results for the gap-center position and the width of the gap are in good agreement with an earlier study by Liew. et al. using a different method but applied to a practically identical system at ϕ = 0.28 volume filling fraction 26 , see also Supplementary Fig. 2. Indeed, for ν 0 2 ½0:47; 0:49, where the DOS is zero, we find that T(L) decays exponentially with L B < a, Fig. 1d. The decay length rises toward the band edges and is smallest around the center frequency ν 0 Gap . This observation is consistent with tunneling of evanescent waves through the whole sample of thickness L 3, 44 . We find that the evanescent regime appears to extend over a slightly broader range of frequencies compared to the bandgap. We will address this point again at the end of this section.
Transmittance in the multiple scattering regime. Standard diffusion theory describes transmission through samples whose thickness L is much larger than the mean free path ℓ and thus T~ℓ/L ≪ 1. Since the size of our simulation box is limited to L ≤ 18a, we also have to include data for slabs with thicknesses on the order of a few ℓ, in particular far from the gap where the transmittance T can be closer to one. To be able to describe the transition from ballistic to diffusive transport as well as to SAL we follow an approach developed by Durian and coworkers. Their theory, which is based on the telegrapher equation, takes into account ballistic transmission, lower order scattering as well as diffusive multiple scattering of light. Their theory accurately predicts T L→0 = 1 for thin samples, while diffusion theory fails in this limit. In the methods section we explain how to consistently merge Durian's approach 52 with the SC theory of localization and we obtain: whereL, D 0 ð Þ and η depend on L; ' ð Þ as well as the localization parameter kl ð Þ. z 0 denotes the extrapolation length in diffusion theory. We note that merging Durian's theory with SC theory is unproblematic. The improvements of the former affect thin slabs L~ℓ while the latter only affects thick slabs L > ξ ≫ ℓ. In essence, Eq. (1) provides a proper interpolation scheme between the two limiting cases. In the absence of SAL, for kl ð Þ ! 1, we recover Durian's results for T L='; z 0 ð Þ withL L, D B /D(0) = 1 and η L=' ð Þ ¼ 0. Then, for L ≫ ℓ, Eq. (1) reduces to the common expression for diffuse transport T L ) '; z 0 ð Þ¼ 1þz 0 2z 0 þL=' . When deriving Eq. (1) we did not distinguish between the scattering and the transport mean free path and only use one ℓ for the mean free path. By considering an extended version of the model, taking into account the scattering anisotropy parameter g, as described in ref. 52 , we have verified that these simplifications do not adversely influence the quality of our fit. Moreover, we neglect specular reflections at the interface of the order of a few percent at most.
In diffusion theory, the extrapolation length z 0 is linked to the angular averaged specular reflectivity R via z 0 ¼ 2 3 1þR 1ÀR . We determine R directly from a fit to the data using a least-squares fitting procedure. Since the transmittance varies by several orders of magnitude and we are interested in the behavior in different regimes, we choose the natural logarithm of the transmittance as the function to fit. For a given frequency we define: as the sum of squares to be minimized. T L i ð Þ is defined according to Eq. (1). The index i runs over the N = 15 data points from L/a = 2.7 to L/a = 18, see Table 1. We note that ℓ sets the optical thickness in terms of the characteristic length a in our system: L/a = L/ℓ × ℓ/a. We first fit Eq. (1) to the numerical data, assuming the absence of localization or kl ð Þ ! 1. Our analysis covers the entire frequency range considered, ν 0 2 ½0:3; 0:6 treating both R and ℓ as adjustable parameters. Figure 2a shows the frequency dependence of S and the values of Rðν 0 Þ we obtain.

Gap
≳ 0:1 we find small S values signaling excellent agreement between diffusing theory and data, as shown in Fig. 2b. The good fit suggests a classical transport regime controlled by ballistic transmission, low order scattering, which eventually evolves to become diffusive for L ≫ ℓ.
In the same frequency range we find R = 0.66 ± 0.05 to be approximately constant, corresponding to z 0 ≃ 3.25 ± 0.5. We repeat the fit keeping R = 0.66 fixed and the goodness of the fit is  Supplementary  Fig. 4). We note that around ν 0 ¼ a=λ Gap $ 0:45 the fit is very poor and the fitted values of ℓ and R become meaningless.
Self-consistent theory of localization. To assess the breakdown of wave diffusion near the band edge, we compare our numerical data to the SC theory of localization. For slabs of finite thickness L we can calculate D(z) for a given value of kl ð Þ 40 . To illustrate the dependence of the diffusion coefficient on z, in Fig. 3 we plot the D as a function of z/L for different values of L/ℓ and kl ð Þ. As expected, for kl ð Þ greater than one the diffusion coefficient shows a weak dependence on both the total size of the system and the position. On the contrary, deep in the localization regime, D(z) decays exponentially away from the boundaries. Exactly at the localization transition for kl ð Þ ¼ 1, the diffusion coefficient is strongly reduced in the center of the slab but the decay remains nonexponential.
We integrate D(z) to obtainL ¼ R L 0 D B =DðzÞ dz and the function η L=' ð Þ. For kl ð Þ values larger than one, SC theory gradually approaches the prediction by diffusion theory and the quality of the fit becomes insensitive to the choice of kl ð Þ. As shown in Fig. 4a, b the S values of both diffusion theory and SC theory become comparable for ν 0 < 0:4 and ν 0 > 0:51, signaling a diffusive transport regime. When approaching the gap from lower (or higher) frequencies, the fit with SC theory, however, leads to substantially smaller S values, indicating localization. From the two-parameter fit we find a lower frequency mobility edge (kl) = 1 at ν 0 c;l ¼ 0:412 (ℓ c,l /a = 0.513) and a higher frequency mobility edge at ν 0 c;h ¼ 0:506 (ℓ c,h /a = 0.242), Fig. 4c, d. Using these anchor values for the mobility edge, (kl) and ℓ are linked via the relation ðklÞ ¼ '=' c ν 0 =ν 0 c . With (kl) ≡ 1 at the mobility edge, we find k ⋅ ℓ = 1.33 for ν 0 c;l ¼ 0:412 and k ⋅ ℓ = 0.77 for ν 0 c;h ¼ 0:506. We note that we use these two separate proportionality constants (1.33, 0.77) for the comparison between theory and numerical data in the higher and lower frequency branch. In all cases, we perform a least-squares fit to ln T according to Eq.
(2). We have considered other tests, such as χ 2 , but these other tests are often based on assumptions that are probably not met in our case. It is for example well known that ln T does not necessarily obey Gaussian statistics in the SAL regime 30 .
Next, we attempt to describe the data with SC theory over the full range of frequencies ν 0 using a single adjustable parameter (kl), as a measure of the distance to the critical point at (kl) = 1. We find excellent agreement between SC theory and the data over the entire range ν 0 ≠½0:45; 0:5, i.e., outside a central frequency interval in or near the full bandgap, Figs. 5 and 6. The observation of such a single-parameter scaling is the key finding of our work. In the regime where (kl) ≲ 1.2 SC theory describes the data substantially better than diffusion theory and we can describe the data across the critical transition from light diffusion to localization. In the low-frequency branch, between the mobility edge ν 0 c;l ¼ 0:412 and the band edge ν 0 ' 0:47, the sample remains localized over a large frequency interval, comparable or larger in width than the bandgap, as shown in Fig. 5. In this regime, the localization parameter drops from (kl) = 1 to about (kl) = 0.85. In the high-frequency branch, the localized regime is nearly degenerate and the sample rapidly enters the gap after crossing the mobility edge from above for ν 0 < ν 0 c;h . For ν 0 < 0:39 and ν 0 > 0:52 the position dependence of D(z) is weak, for our system sizes, and the predictions by SC theory and by diffusion theory are indistinguishable.

Discussion
In summary, we could show that hyperuniform 3D silicon networks display different characteristic transport regimes for electromagnetic vector waves. Deep in the gap region, the transmittance decays exponentially, indicating tunneling through the entire sample. Outside the gap region, we observe a critical transition from classical diffusion to wave localization controlled by a single parameter (kl). We have shown that in this regime, our numerical data can be described quantitatively by assuming a position-dependent diffusion coefficient D(z) derived from the SC theory of localization. Finding such an agreement is generally understood as evidence for SAL of light. Moreover, our detailed results about the transition to SAL in realistic digital representations of dielectric networks can provide valuable guidance for future experimental attempts to probe light localization.
Finally, we would like to add a remark concerning the breakdown of the description of wave transport by SC theory and the opening of the gap. SC theory is a heuristic approach that postulates a position-dependent diffusion coefficient D(z) but does not provide a microscopic explanation (which, under certain conditions has been added later 38 ). The concept of applying a diffusion equation to describe transport over certain distances ξ > ℓ breaks down as ξ/ℓ → 1. It it thus unclear whether evanescent decay, Fig. 1d, and the correspondingly poor fit with SC theory, within Δν 0 ' 0:01 À 0:02 next to the band edge, is due to the breakdown of SC theory or other emerging transport phenomena in the vicinity of the gap, as suggested earlier in ref. 30 . Moreover, we find it interesting to speculate whether the SC theory could be generalized to provide a unified theory for wave transport in amorphous photonic materials encompassing the transition between diffusive, localized, and bandgap regimes.

Methods
FDTD simulations. FDTD simulations were performed using the MEEP software 49 and were run on a computer cluster. Throughout this work, the Poynting vector is recorded on a monitor situated behind the structure. Transmittance is defined as the ratio of the transmitted power to the incoming one and is calculated by dividing the transmitted power (integral of the Poynting vector over the monitor) by the power transmitted in a reference run (empty simulation box). The network structures were generated using a custom-made code (MATLAB and Statistics Toolbox, The MathWorks, Inc., Massachusetts, United States) based on the full 10,000 particle seed pattern taken from 48 . Equivalently, the sicipy.spatial opensource library of Python can be used for this purpose. Using a clean cut we obtain slabs of different thickness L ≤ 18a, which were then imported to MEEP (see Table 1 for the exact values of thickness). All units were set to μm and the sphere diameter was set to a = 5/3 μm. We apply periodic boundary conditions in x-and y-directions. Since we do not know the precise size of the simulation box used in ref. 48 , the periodic boundary conditions in MEEP will not exactly match the original periodic boundary conditions employed when generating the pattern. From our own band structure calculations we find that this can give rise to a few  defect states in the gap. We do not expect these defect states to contribute to transport in the diffuse or SAL regime, where the DOS is high, but they can increase the transmittance in the gap for L ≫ L B due to tunneling between defect states. We believe this is the main reason why TðLÞ> expðÀL=L B Þ for L ≫ L B in the gap regime, see e.g., Fig. 6d. Since the gap regime is not in the focus of our study, we have not explored this in more detail but plan to address this in future work, using new seed pattern properly matched to the periodic boundary conditions of the MEEP-simulation box.
The network was illuminated with a broadband pulse of linearly polarized light with electric field vector parallel to one of the sides of the simulation box. The pulse bandwidth was sufficient to cover vacuum wavelengths between 1 and 7 μm corresponding to reduced frequencies ν 0 2 ½0:24; 1:67. Perfectly matched layers (PML) were fitted at both ends of the simulation box and periodic boundary conditions were applied perpendicular to the wave propagation direction. The PML's absorb all transmitted and reflected waves (regardless of incidence direction) and prevent them from re-entering the simulation box. The PML thickness was 7 microns, which ensured that all wavelengths shorter than this value were suppressed. The spatial resolution was equal to 20 pixels per μm. Convergence tests were performed to check the robustness of the simulation. It was verified that increasing the spatial resolution by a factor of two did not considerably influence the transmittance curves. Also the simulation time was selected in a way to yield robust results. By placing an additional monitor between the source and the network, we checked for flux conservation over the entire frequency interval of interest.
Numerical implementation of the SC theory of localization. In the selfconsistent theory of localization, the standard diffusion equation is replaced by an equation where the diffusion coefficient is nonlocal both in the space and time domain. The renormalization of the diffusion coefficient accounts for the different return probability when interference effects are considered in the multiple scattering regime 38 . We are considering a slab geometry and continuous wave illumination at a given carrier wave frequency. The simplified geometry, invariant in the plane parallel to the slab, leads to a set of SC equations for the diffusion constant dependence and the diffusion equation Green's function g z; z 0 ð Þ in the direction z perpendicular to the slab surfaces (lying between z = 0 and z = L). Here we work in reduced units, where all lengths are scaled by the mean free path ℓ. The diffusion coefficient D z ð Þ is normalized by the Botzmann diffusion coefficient D B . Taking the Fourier transform in the x, y plane, parallel to the slab boundaries, we reach at the diffusion equation: This result together with the self-consistent equation for the diffusion coefficient determines the transport properties of our system for all values of the parameters. In the above equation the cut off is given by q max ¼ 1=3ðklÞ 2 c . The latter depends on the chosen value (kl) c with (kl) = (kl) c at the mobility edge. We solve the set of self-consistent equations by recursively solving Eq. (3) in a first step for all positions of the source z 0 and transversal wavenumbers q. The solution of the Green function is then plugged into Eq. (4) to correct for the zdependent diffusion coefficient D z ð Þ which is inserted back to Eq. (3) until convergence is reached. We choose a discretization scheme where all the considered positions z and wavenumbers squared q 2 are evenly spaced. Depending on the values of the parameters (kl), (kl) c , and total length L, we need to take a step Δz ≡ h ranging from 0.02 to 0.2 and between 300 and 600 steps in q 2 to achieve a relative precision in D z ð Þ of the order or 10 −4 in 5-50 recursion steps. Taking the second order finite differences approximation for the derivatives in Eq. (3) leads to a tridiagonal system of equations which has to be solved for each value of the wavenumber q and position of the source. Obviously, changing the position of the source amounts to changing the independent term of the system of equations and hence all equations for a given value of q can be solved at once through the inverse of the corresponding tridiagonal matrix. We choose the Lapack function DGTSV to get the inverse since it is simple to use and universally accessible while efficient enough for our purposes 53 . Specifically, if we take a discretization , the diagonal terms of the system of equations read the subdiagonal terms are Analogously, the superdiagonal terms are Since the sources are located in the interior of the slab, the points at which the source is located are z 0 i¼2; ÁÁÁ; nÀ1 ¼ h i À 1 ð Þ. The independent term vector T i for the source at z i is T i ð Þ j ¼ hδ i;j . Once the numerical solution is obtained for all the source positions, the value of g z i ; z i ð Þis available for i = 2, ⋯, n − 1. The g z 1 ; z 1 ð Þand g z n ; z n ð Þare obtained by second order accuracy extrapolation of the solutions g z i ; z i ð Þto the boundary. The full solution is then used to calculate the integral in Eq. (4) by 3/8 Simpson's rule. The updated value of the diffusion coefficient D z ð Þ is used in the next step of recursion. The algorithm stops when the logarithm of an update of D z ð Þ differs by <10 −4 from the previous value at all the points in the discretization. The seed D z ð Þ used to start the recursions is set to 1. In all cases (kl) c = 1.
Transmission through a slab. In the following, we describe how to consistently merge SC theory with the theory by Lemieux et al. 52 . The contribution that is affected by SC theory is the total diffusive transmittance T d (see also 54 ). We compute it considering that the sources of diffuse intensity are continuously distributed across the slab with an intensity proportional to the ballistic intensity I b ¼ exp Àz ð Þ. To simplify the notation, in this section, we work in reduced units where all lengths are normalized to ℓ and D(z) is normalized to D B . The diffusive transmittance, neglecting boundary reflectivity of a few percent, is hence where T SCT z ð Þ is the diffuse transmittance for a source at z ¼ z 0 (see Supplementary Material, Eqs. (S5)-(S17), for details), explicitly we have withL The first term (z 0 in the numerator) can be explicitly integrated: The second term in Eq. (9) is more involved since D z ð Þ is not known a priori, the integral to be solved is that can be formally integrated by parts using u z ð Þ R z 0 1 D x ð Þ dx and dv z ð Þ e Àz dz. Hence The second term in the rhs of Eq. (13) can again be formally integrated by parts using u z ð Þ 1=DðzÞ, dv z ð Þ e Àz dz to give Collecting results, we have Finally, adding the ballistic transmittance e −L , we obtain Eq. (1) for the total transmittance: It is worth noticing at this point that in the limit of standard diffusion theory, D z ð Þ ¼ 1 (i.e. the diffusion coefficients coincides with the standard D B = cℓ/3) and η L ð Þ ¼ 0, since LN z ð Þ ¼ 0).

Data availability
The data of T(L/a) shown in the paper, either obtained from FDTD calculations or from the self-consistent theory of localization (SC theory) are provided in the online repository (https://doi.org/10.5281/zenodo.3968424). These and all other data sets used in the study can be generated from the codes uploaded to the repository, or can be obtained from us upon reasonable request.

Code availability
The codes used to produce the results of this study are included in the repository https:// doi.org/10.5281/zenodo.3968424 and described in the main text or the supplementary material. With respect to third party codes, such as the open-source codes MPB and MEEP, we refer to the original publications, see refs. [49][50][51] . Links to the original third party sources are provided in the README.txt files in the repository.