\newcites

suppSupplementary References

Anderson localization of electromagnetic waves in three dimensions

Alexey Yamilov^1⋆, Sergey E. Skipetrov², Tyler W. Hughes³, Momchil Minkov³,
Zongfu Yu^3,4,†, Hui Cao^5∗

¹Physics Department, Missouri University of Science & Technology, Rolla, Missouri 65409
²Univ. Grenoble Alpes, CNRS, LPMMC, 38000 Grenoble, France
³Flexcompute Inc, 130 Trapelo Road, Belmont, MA, 02478
⁴Department of Electrical & Computer Engineering, University of Wisconsin, Madison, WI, 53705
⁵Department of Applied Physics, Yale University, New Haven, Connecticut 06520
^⋆ yamilov@mst.edu
^† zongfu@flexcompute.com
^∗ hui.cao@yale.edu

Anderson localization marks a halt of diffusive wave propagation in disordered systems. Despite extensive studies over the past 40 years, Anderson localization of light in three dimensions has remained elusive, leading to the question of its very existence. Recent orders-of-magnitude speed-up of finite-difference time-domain calculations allows us to conduct brute-force numerical simulations of light transport in fully disordered 3D systems with unprecedented dimension and refractive index contrast. We demonstrate three-dimensional localization of vector electromagnetic waves in random packings of metallic spheres, in sharp contrast to the absence of localization for dielectric spheres with a refractive index contrast up to 10. Our work opens a wide range of avenues in both fundamental research related to Anderson localization and potential applications using 3D localized light.

Anderson localization (AL) (?) is an emergent phenomenon for both quantum and classical waves including electron (?, ?, ?), cold-atom (?, ?), electromagnetic (EM) (?, ?, ?, ?, ?), acoustic (?, ?), water (?), seismic (?), and gravity (?) waves. Unlike in one or two dimensions, AL in three dimensions (3D) requires strong disorder (?, ?, ?, ?). A mobility edge separating the diffuse transport regime from AL can be estimated from the Ioffe-Regel criterion $k_{\text{eff}}\,\ell_{s}\sim 1$ , where $k_{\text{eff}}$ is the effective wavenumber in the medium and $\ell_{s}$ is the scattering mean free path (?). This criterion suggests two avenues to achieving localization: reduction of $k_{\text{eff}}$ or $\ell_{s}$ . For EM waves, the former is realized by introducing partial order or spatial correlation in scatterer positions (?, ?). In comparison, reaching localization of light in fully random photonic media by increasing scattering strength (decreasing $\ell_{s}$ ) turns out to be much more challenging (?, ?). Despite successful experiments in low-dimensional systems (?, ?, ?), 3D localization remained stubbornly elusive (?), which triggered theoretical (?, ?) and experimental (?) studies of the mechanisms that impede it.

Anderson himself originally proposed “a system composed essentially of random waveguides near cut-off and random resonators, such as might be realized by a random packing of metallic balls of the right size” as “the ideal system” for localization of EM radiation (?). In practice, the absorption of metals obscures localization (?, ?), and the experimental focus shifted to dielectric materials with low loss and high refractive index (?, ?, ?, ?, ?). However, even for dielectric systems, experimental artifacts due to residual absorption and inelastic scattering mar the data (?, ?). Numerically, these artifacts can be excluded, but 3D random systems of sufficiently large dimension and high refractive index contrast could not be simulated due to an extraordinarily long computational time required (?, ?).

Recent implementation of the finite-difference time-domain (FDTD) algorithm on emerging computing hardware has brought orders of magnitude speed-up of numerical calculation (?, ?). Using this highly-efficient hardware-optimized version of the FDTD method, we solve the Maxwell equations by brute force in 3D. It enables us to simulate sufficiently large systems and high refractive index contrast to address the following questions: can 3D AL of EM waves be achieved in fully random systems of dielectric scatterers? If not, can it occur in any other systems without aid of spatial correlations?

Answering these long-standing questions not only addresses the fundamental aspects of wave transport and localization across multiple disciplines, but also opens new research avenues and applications. For example, in topological photonics (?), the interplay between disorder and topological phenomena may be explored beyond the limit of weak disorder in low-dimensional systems (?). Also in cavity quantum electrodynamics with Anderson-localized modes (?), going 3D will avoid the out-of-plane loss inherent for 2D systems and cover the full angular range of propagation directions (?). In addition to fundamental studies, disorder and scattering has been harnessed for various photonic device applications, but mostly with diffuse waves (?). Anderson-localized modes can be used for 3D energy confinement to enhance optical nonlinearities, light-matter interactions, and control random lasing as well as targeted energy deposition.

We first consider EM wave propagation through a 3D slab of randomly packed lossless dielectric spheres of radius $r=100$ nm and refractive index $n=3.5$ . This corresponds to the highest index contrast achieved experimentally in optical range with porous GaP around the wavelength $\lambda=650$ nm in the vicinity of the first Mie resonance of an isolated sphere, see Supplementary Information (SI). To avoid spatial correlations, the sphere positions are chosen completely randomly, leading to spatial overlap where index is capped at the same value of $n$ . We compute the spatial correlation function of such structure, which reveals that the correlation vanishes beyond the particle diameter (SI). To avoid light reflection at the interfaces of the slab, we surround it by a uniform medium with the refractive index equal to the effective index of the slab, $n_{\text{eff}}=[(1-f)+fn^{2}]^{1/2}$ , for a given dielectric volume filling fraction $f$ (Fig. 1a). As described in SI, for each wavelength, we compute the scattering mean free path $\ell_{s}$ directly from the rate of attenuation of co-polarized field with depth. This, together with the effective wavenumber $k_{\text{eff}}=n_{\text{eff}}\,(2\pi/\lambda)$ , gives the Ioffe-Regel parameter plotted in Fig. 1c. It features a minimum around $\lambda=650$ nm, and the smallest value of $k_{\text{eff}}\,\ell_{s}\simeq 0.9$ is reached at $f=38\%$ . We also compute the transport mean free path $\ell_{t}$ from continuous wave (CW) transmittance of an optically thick slab with thickness $L\gg\ell_{t}$ (SI). In Fig. 1d, $k_{\text{eff}}\,\ell_{t}$ also exhibits a dip in the same wavelength range as $k_{\text{eff}}\,\ell_{s}$ , however, the smallest $k_{\text{eff}}\,\ell_{t}$ is found at lower $f$ of $18$ –29%, as the dependent scattering sets in at higher $f$ . In search for AL in this wavelength range, we numerically simulate propagation of a narrowband Gaussian pulse centered at $\lambda_{0}=650$ nm with plane wavefront, and compute the transmittance through the slab $T(t)$ as a function of arrival time $t$ . The diffusive propagation time $\tau_{D}$ approximately corresponds to the arrival time of the peak in Fig. 1e. At $t\gg\tau_{D}$ , the decay of the transmitted flux is exponential over $12$ orders of magnitude, as expected for purely diffusive systems (?). Rate of this exponential decay is equal to $1/\tau_{D}$ , which is directly related (?) to the smallest diffusion coefficient within the spectral range of the excitation pulse (SI). In Fig. 1f, the dependence of this diffusion coefficient $D$ on the dielectric filling fraction $f$ exhibits a minimum at $f\sim 30\%$ . Inset in Fig. 1e shows that the further increase of the slab thickness does not lead to any deviation from diffusive transport. Furthermore, the diffusive behavior persists in the numerical simulation with increased spatio-temporal resolution (SI). At $t\gg\tau_{D}$ , the spatial intensity distribution inside the system features a depth profile (averaged over cross-section) equal to that of the first eigenmode of the diffusion equation (Fig. 1b). We therefore rule out a possibility of AL in uncorrelated ensembles of dielectric spheres with $n=3.5$ .

At microwave frequencies, the refractive index may be even higher than $n=3.5$ . We therefore, perform numerical simulation of a 3D slab of dielectric spheres with $n=10$ . The main results are summarized here, and details are presented in SI. A large scattering cross section $\sigma_{s}(\lambda)$ of a single sphere near the first Mie resonance leads to strong dependent scattering already at small filling fractions. We find the Ioffe-Regel parameter $k_{\text{eff}}\,\ell_{s}\gtrsim 1$ despite the very high refractive index contrast. This is attributed to dependent scattering that becomes significant even at relatively low dielectric filling fraction $f$ . The numerically calculated $T(t)$ for $L/\ell_{t}\gg 1$ does not exhibit any deviation from diffusive transport: at $t\gg\tau_{D}$ , the decay of transmittance is still exponential over $\sim 10$ orders of magnitude. In addition, scaling of CW transmittance with the inverse slab thickness $1/L$ remains linear for all $f$ , as expected for diffusion (SI). We therefore conclude that AL does not occur in random ensembles of dielectric spheres, thus closing the debate about the possibility of light localization in a white paint (?, ?).

Previous studies (?, ?) suggest that absence of AL for EM waves may be due to longitudinal waves that exist in a heterogeneous dielectric medium, where the transversality condition $\nabla\cdot\mathbf{E}(\mathbf{r})=0$ does not follow from the Gauss’s law $\nabla\cdot[\epsilon(\mathbf{r})\mathbf{E}(\mathbf{r})]=0$ , because of the position dependence of $\epsilon(\mathbf{r})$ . Here, we propose to suppress the contribution of longitudinal waves to optical transport and realize AL of EM waves by using perfectly conducting spheres as scatterers. The Poynting vector is parallel to the surface of a perfect electric conductor (PEC) (?) and EM energy flows around a PEC particle without coupling to longitudinal surface modes. The volume of PEC spheres is simply excluded from the free space and becomes unavailable for light. Thus, at high PEC volume fraction, light propagates in a random network of irregular air cavities and waveguides formed by the overlapping PEC spheres, akin to the original proposal of Anderson (?).

Similarly to the dielectric systems above, we simulate a 3D slab composed of randomly packed, overlapping PEC spheres of radius $r=50$ nm in air. Figure 2 shows the results of simulating an optical pulse propagating through 10 $\mu$ m $\times$ 10 $\mu$ m $\times$ 3.3 $\mu$ m slabs of various PEC fractions. $T(t)$ displays non-exponential tails at high fractions $f=41\%,\ 48\%$ in Fig. 2a. From the decay rate obtained via a sliding-window fit, we extract a time-dependent diffusion coefficient $D(t)$ , c.f. Fig. 2b, which shows a power-law decay with time, as predicted by the self-consistent theory of localization (?). The non-exponential decay of $T(t)$ and the time-dependence of $D$ are the signatures of AL (?, ?). In contrast, at lower PEC fractions $f=8\%,\ 15\%$ , $D$ remains constant in time. Figure 2c reveals a transition from time-invariant $D$ to time-dependent $D(t)$ around $f=33\%$ where $D(t)$ starts deviating from a constant. Using Fourier transform, we compute the spectrally resolved transmittance $T(\lambda)$ . Figures 2d,e contrast $T(\lambda)$ of diffusive and localized systems. The former features smooth, gradual variations with $\lambda$ due to broad overlapping resonances, whereas the latter exhibits strong resonant structures consistent with the average mode spacing exceeding the linewidth of individual modes, in accordance with Thouless criterion of localization, as the spectral narrowing of modes is intimately related to their spatial confinement (?, ?, ?). Color maps in Figs. 2d,e show spatial intensity distributions inside the systems, $\langle I(x,y_{0},z;\lambda)\rangle_{x}$ , averaged over $x$ for a cross-section $y=y_{0}$ . These two-dimensional maps contrast slow variation with $z$ and $\lambda$ in the diffusive system (panel d) to the sharp features due to spatially-confined modes in the localized system (panel e). Furthermore, there exist ‘necklace’ states with multiple spatially separated intensity maxima, originally predicted for electrons in metals (?).

Insight into the mechanism behind AL in the random ensemble of PEC spheres can be gained from the wavelength dependence of the Ioffe-Regel parameter $k\ell_{s}$ . We compute it using a procedure similar to that in dielectrics (SI). Even at the volume fraction of $f=8\%$ , $\ell_{s}$ is well below the prediction of independent scattering approximation (ISA), due to scattering resonances formed by two or more neighboring PEC spheres (SI). As shown in Fig. 3a, $\ell_{s}$ becomes essentially independent of wavelength in the range of size parameter $kr$ of PEC spheres. Consequently, the Ioffe-Regel parameter acquires $1/\lambda$ dependence as seen in Fig. 3b. It drops below the value of unity within the excitation pulse bandwidth $\lambda\simeq 650$ nm $\pm$ 45 nm for $f$ between $25\%$ and $33\%$ , in agreement with Fig. 2. We further conduct a finite-size scaling study, after computing the dependence of CW transmittance $T$ on the slab thickness $L$ (SI). Figure 3c shows the logarithmic derivative $d\log(T)/d\log(L)$ as a function of $k\ell_{s}$ . In the diffusive regime, Ohm’s law $T\propto 1/L$ is expected, leading to a scaling power of $-1$ , as indeed observed for $k\ell_{s}>1$ . Around $k\ell_{s}\sim 1$ we see a departure from $1/L$ scaling of transmittance. Note that the scaling theory of localization (?) predicts a single-parameter scaling with the dimensionless conductance $g$ but not with $k\,\ell_{s}$ . By estimating the number of transverse modes as $N=2\pi(L/\lambda)^{2}(1-f)^{2/3}$ for $L\times L$ area of the slab, we compute $g=TN$ (?) and $\beta(g)\equiv d\log(g)/d\log(L)$ . Figure 3d shows a good agreement between the numerical data and the model function $\beta(g)=2-(1+g)\log(1+g^{-1})$ (?). In diffusive regime $g>1$ , $\beta(g)\rightarrow 1$ ; whereas in the localized regime $g<1$ , $\beta(g)\propto\log(g)$ . The latter is a manifestation of the negative exponential scaling of $g$ with $L$ in the regime of Anderson localization.

To obtain the ultimate confirmation of AL of light in PEC composites, we simulate dynamics of the transverse spreading of a tightly focused pulse—a measurement that has been widely adopted in localization experiments (?, ?, ?). A pulse centered at $\lambda=650$ nm with a bandwidth of 90 nm is focused to a small spot of area $\simeq 0.5$ $\mu\rm{m}^{2}$ at the front surface of a wide 3D slab of dimensions 33 $\mu$ m $\times$ 33 $\mu$ m $\times$ 3.3 $\mu$ m, c.f. Fig. 4a. We compute the transverse extent of intensity distribution $I(x,y,z=L;t)$ at the back surface of the slab. For a diffusive PEC slab with $f=15\%$ , we observe a rapid transverse spreading of light with time in Fig. 4b, which approaches the lateral boundary of the slab within $\sim$ 2 ps. In sharp contrast, in the localized system in Fig. 4c ( $f=48$ %), the transmitted intensity profile remains transversely confined even after 20 ps. This time corresponds to a free space propagation of 6 mm, which is $\sim 2000$ times longer than the actual thickness of the slab. Figure 4d quantifies this time evolution with the output beam diameter $d(t)=2[\text{PR}(t)/\pi]^{1/2}$ , where $\text{PR}(t)=[\int\int I(x,y,L;t)dxdy]^{2}/\int\int I(x,y,L;t)^{2}dxdy$ is the intensity participation ratio. For a diffusive slab, $d(t)\propto t^{1/2}$ , while in the localized regime, $d(t)$ saturates at a value of the order of slab thickness $L$ . Such an arrest of the transverse spreading in the localized PEC systems persists with increased spatio-temporal resolution of the numerical simulation (SI). Further evidence of AL include non-linear depth profile and strong non-Gaussian fluctuations of intensity inside the system (SI). We also confirm our results by repeating calculations for 3D slabs of PEC spheres with larger radius $r=100$ nm, and obtaining similar scaling behavior (SI), as in Figs. 3cd.

The striking difference between light propagation in dense random ensembles of dielectric and PEC spheres cannot be accounted for by the Ioffe-Regel parameter, as both reach $k_{\text{eff}}\,\ell_{s}\sim 1$ for similar values of the size parameter $kr$ (Figs. 1c, 3b). AL in 3D PEC composites with uncorrelated disorder reveals a localization mechanism that is unique for PEC. In contrast to a dielectric system where light propagates everywhere (both inside and outside the scatterers), the propagation is restricted to the voids between scatterers in the PEC system. AL takes place when the typical size of voids becomes smaller than the wavelength, and light can no longer ‘squeeze’ through them. This qualitative model correctly predicts both AL in the long-wavelength regime and the increase of the critical volume fraction $f$ for localization with the scatterer size $r$ (SI).

Finally, we test AL in real-metal aggregates. In the microwave spectral region, the skin depth of crystalline metals like silver, aluminum and copper is several orders of magnitude shorter than the wavelength $\lambda$ and the scatterer size $r$ in the regime of $kr\sim 1$ . Since the microwave barely penetrates into the metallic scatterers, our simulation results are almost identical to those for PEC. To account for the imperfections due to polycrystallinity, surface defects, oxide layers etc., we lower the metal conductivity to match the experimentally measured absorption rate in aggregates of aluminum spheres (?). Simulations unambiguously show the arrest of transverse spreading of a focused pulse (Fig. 4e), demonstrating AL in 3D random aggregates of aluminum spheres. Additional evidence of AL is presented in SI. Moreover, even at optical frequencies, where realistic metals deviate notably from PEC, the arrest of transverse spreading persists in 3D silver nanocomposites, as shown in SI. Possible light localization in 3D nanoporous metals will have a profound impact on their applications in photo-catalysis, optical sensing, energy conversion and storage.

In summary, our large-scale microscopic simulations of EM wave propagation in 3D uncorrelated random ensembles of spheres show no signs of Anderson localization for dielectric spheres with refractive indices $n=3.5$ –10. This may explain multiple failed attempts of experimental observation of AL of light in 3D dielectric systems over the last three decades (?, ?, ?, ?, ?). At the same time, we report the first numerical evidence of EM wave localization in random ensembles of PEC spheres over a broad spectral range. Localization is confirmed by eight criteria: the Ioffe-Regel criterion, the Thouless criterion, non-exponential decay of transmittance under pulsed excitation, vanishing of the diffusion coefficient, existence of spatially localized states, scaling of conductance, arrest of the transverse spreading of a narrow beam, and enhanced non-Gaussian fluctuations of intensity. Our study calls for renewed experimental efforts to be directed at low-loss metallic random systems (?). In the SI, we propose a realistic microwave experiment that avoids experimental pitfalls (?) and provides a tell-tale sign of Anderson localization.

References

1. Anderson, P. W. Absence of diffusion in certain random lattices. Phys. Rev. 109, 1492–1505 (1958).
2. Mott, N. Electrons in disordered structures. Adv. Phys. 16, 49–144 (1967).
3. Kramer, B. & MacKinnon, A. Localization: theory and experiment. Rep. Prog. Phys. 56, 1469–1564 (1993).
4. Imada, M., Fujimori, A. & Tokura, Y. Metal-insulator transitions. Rev. Mod. Phys. 70, 1039–1263 (1998).
5. Billy, J. et al. Direct observation of Anderson localization of matter waves in a controlled disorder. Nature 453, 891–894 (2008).
6. Jendrzejewski, F. et al. Three-dimensional localization of ultracold atoms in an optical disordered potential. Nat. Phys. 8, 398–403 (2012).
7. John, S. Electromagnetic absorption in a disordered medium near a photon mobility edge. Phys. Rev. Lett. 53, 2169–2172 (1984).
8. Anderson, P. W. The question of classical localization. A theory of white paint? Phil. Mag. B 52, 505–509 (1985).
9. Chabanov, A. A., Stoytchev, M. & Genack, A. Z. Statistical signatures of photon localization. Nature 404, 850–853 (2000).
10. Schwartz, T., Bartal, G., Fishman, S. & Segev, M. Transport and Anderson localization in disordered two-dimensional photonic lattices. Nature 446, 52–55 (2007).
11. Segev, M., Silberberg, Y. & Christodoulides, D. N. Anderson localization of light. Nat. Photonics 7, 197–204 (2013).
12. Kirkpatrick, T. R. Localization of acoustic waves. Phys. Rev. B 31, 5746–5755 (1985).
13. Hu, H., Strybulevych, A., Page, J. H., Skipetrov, S. E. & van Tiggelen, B. A. Localization of ultrasound in a three-dimensional elastic network. Nat. Phys. 4, 945–948 (2008).
14. Guazzelli, E., Guyon, E. & Souillard, B. On the localization of shallow water waves by random bottom. J. Phys. Lett. 44, 837–841 (1983).
15. Sheng, P., White, B., Zhang, Z. Q. & Papanicolaou, G. Wave localization and multiple scattering in randomly-layered media. In Sheng, P. (ed.) Scattering and Localization of classical waves in random media, Directions in Condensed Matter Physics, 563–619 (World Scientific Publishing, 1990).
16. Rothstein, I. Z. Gravitational Anderson localization. Phys. Rev. Lett. 110, 011601 (2013).
17. John, S. Localization of light. Phys. Today 44, 32–40 (1991).
18. Sheng, P. Introduction to wave scattering, localization and mesoscopic phenomena (Springer, 2006).
19. Lagendijk, A., van Tiggelen, B. & Wiersma, D. S. Fifty years of Anderson localization. Phys. Today 62, 24–29 (2009).
20. Ioffe, A. F. & Regel, A. R. Non-crystalline, amorphous, and liquid electronic semiconductors. Prog. Semicond. 4, 237–291 (1960).
21. Haberko, J., Froufe-Perez, L. S. & Scheffold, F. Transition from light diffusion to localization in three-dimensional amorphous dielectric networks near the band edge. Nat. Comm. 11, 4867 (2020).
22. van der Beek, T., Barthelemy, P., Johnson, P. M., Wiersma, D. S. & Lagendijk, A. Light transport through disordered layers of dense gallium arsenide submicron particles. Phys. Rev. B 85, 115401 (2012).
23. Sperling, T. et al. Can 3D light localization be reached in ‘white paint’? New J. Phys. 18, 013039 (2016).
24. Lahini, Y. et al. Anderson localization and nonlinearity in one-dimensional disordered photonic lattices. Phys. Rev. Lett. 100, 013906 (2008).
25. Skipetrov, S. E. & Page, J. H. Red light for Anderson localization. New J. Physics 18, 021001 (2016).
26. Skipetrov, S. E. & Sokolov, I. M. Absence of Anderson localization of light in a random ensemble of point scatterers. Phys. Rev. Lett. 112, 023905 (2014).
27. van Tiggelen, B. A. & Skipetrov, S. E. Longitudinal modes in diffusion and localization of light. Phys. Rev. B 103, 174204 (2021).
28. Cobus, L. A., Maret, G. & Aubry, A. Crossover from renormalized to conventional diffusion near the three-dimensional anderson localization transition for light. Phys. Rev. B 106, 014208 (2022).
29. Genack, A. Z. & Garcia, N. Observation of photon localization in a three-dimensional disordered system. Phys. Rev. Lett. 66, 2064–2067 (1991).
30. Watson Jr, G., Fleury, P. & McCall, S. Searching for photon localization in the time domain. Phys. Rev. Lett. 58, 945 (1987).
31. Wiersma, D. S., Bartolini, P., Lagendijk, A. & Righini, R. Localization of light in a disordered medium. Nature 390, 671–673 (1997).
32. Störzer, M., Gross, P., Aegerter, C. M. & Maret, G. Observation of the critical regime near Anderson localization of light. Phys. Rev. Lett. 96, 063904 (2006).
33. Sperling, T., Bührer, W., Aegerter, C. M. & Maret, G. Direct determination of the transition to localization of light in three dimensions. Nat. Photonics 7, 48–52 (2013).
34. Gentilini, S., Fratalocchi, A., Angelani, L., Ruocco, G. & Conti, C. Ultrashort pulse propagation and the Anderson localization. Opt. Lett. 34, 130–132 (2009).
35. Pattelli, L., Egel, A., Lemmer, U. & Wiersma, D. S. Role of packing density and spatial correlations in strongly scattering 3D systems. Optica 5, 1037–1045 (2018).
36. Flexcompute inc. https://flexcompute.com, Accessed: 2021-12-31.
37. Hughes, T. W., Minkov, M., Liu, V., Yu, Z. & Fan, S. A perspective on the pathway toward full wave simulation of large area metalenses. Appl. Phys. Lett. 119, 150502 (2021).
38. Lu, L., Joannopoulos, J. D. & Soljacic, M. Topological photonics. Nat. Photonics 8, 821–829 (2014).
39. Stützer, S. et al. Photonic topological Anderson insulators. Nature 560, 461–465 (2018).
40. Sapienza, L. et al. Cavity quantum electrodynamics with Anderson-localized modes. Science 327, 1352–1355 (2010).
41. Wiersma, D. S. Random quantum networks. Science 327, 1333–1334 (2010).
42. Cao, H. & Eliezer, Y. Harnessing disorder for photonic device applications. Appl. Phys. Rev. (2022). To be published.
43. Akkermans, E. & Montambaux, G. Mesoscopic Physics of Electrons and Photons (Cambridge University Press, Cambridge, UK, 2007).
44. van de Hulst, H. C. Light scattering by small particles (Dover, New York, 1981).
45. Skipetrov, S. E. & van Tiggelen, B. A. Dynamics of Anderson localization in open 3D media. Phys. Rev. Lett. 96, 043902 (2006).
46. Thouless, D. J. Electrons in disordered systems and the theory of localization. Phys. Rep. 13, 93–142 (1974).
47. Pendry, J. B. Quasi-extended electron states in strongly disordered systems. J. Phys. C 20, 733–742 (1987).
48. Abrahams, E., Anderson, P. W., Licciardello, D. C. & Ramakrishnan, T. V. Scaling theory of localization: Absence of quantum diffusion in two dimensions. Phys. Rev. Lett. 42, 673–676 (1979).
49. Müller, C. A. & Delande, D. Disorder and interference: localization phenomena. In Ultracold Gases and Quantum Information: Lecture Notes of the Les Houches Summer School in Singapore, chap. 9 (Oxford University Press, 2011).
50. Cherroret, N., Skipetrov, S. E. & van Tiggelen, B. A. Transverse confinement of waves in three-dimensional random media. Phys. Rev. E 82, 056603 (2010).

Refer to caption — Figure 1: Absence of non-diffusive transport in random dielectric systems with index contrast of 3.5. a, 3D slab filled with dielectric spheres at random uncorrelated positions (radius $r=100$ nm, refractive index $n=3.5$ ) in air. The slab cross-section is 10 $\mu$ m $\times$ 10 $\mu$ m = 100 $\mu$ m² and thickness is $L=3.3$ $\mu$ m. b, 3D distribution of light intensity inside the slab (dielectric filling fraction $f=29$ %, $L/\ell_{t}=33$ ) at long delay time after a short pulse of plane wave front is incident on the front surface. Red curve with shading shows the average depth profile. c, Spectral dependence of the Ioffe-Regel parameter $k_{\text{eff}}\,\ell_{s}$ for different volume filling fractions of dielectric spheres, showing enhancement of scattering around single-sphere Mie resonances. The horizontal dashed line marks the Ioffe-Regel criterion $k_{\text{eff}}\,\ell_{s}=1$ for 3D localization. d, Transport mean free path $\ell_{t}$ (in units of $1/k_{\text{eff}}$ ) as a function of wavelength, revealing a saturation by dependent scattering at high dielectric filling fractions. The vertical dashed lines mark spectral width (33 nm) of the excitation pulse in b and e. e, Transmittance of the 3D slab for a pulsed excitation, showing exponential decay in time for all dielectric filling fractions, in agreement with diffusive transport. Inset shows persistence of diffusion when $L/\ell_{t}$ is increased from 33 to 60 for $f=38\%$ , c.f. green line. f, Dependence of the minimum diffusion coefficient within the pulse bandwidth on the dielectric filling fraction $f$ , exhibiting a minimum value of 3.6 m²/s at $f\simeq 30\%$ .

Acknowledgments

This work is supported by the National Science Foundation under Grant Nos. DMR-1905465, DMR-1905442, and the Office of Naval Research (ONR) under Grant No. N00014-20-1-2197. We thank Shanhui Fan and Bart van Tiggelen for enlightening discussions.

Conflict of interest statement

TWH, MM, ZY have financial interest in Flexcompute Inc., which develops the software Tidy3D used in this work.

Author contributions

A.Y. performed numerical simulations, analyzed the data and compiled all results. S.E.S. conducted theoretical study and guided data interpretation. T.W.H. and M.M. implemented the hardware-accelerated FDTD method and aided in the setup of the numerical simulations. Z.Y. and H.C. initiated this project and supervised the research. A.Y. wrote the first draft, S.E.S. and H.C. revised the content and scope, T.W.H. and M.M. and Z. Y. edited the manuscript. All co-authors discussed and approved the content.

Supplementary Information:
Anderson localization of electromagnetic waves in three dimensions

Alexey Yamilov^1⋆, Sergey E. Skipetrov², Tyler W. Hughes³, Momchil Minkov³,
Zongfu Yu^3,4,†, Hui Cao^5∗

¹Physics Department, Missouri University of Science & Technology,Rolla, Missouri 65409
²Univ. Grenoble Alpes, CNRS, LPMMC, 38000 Grenoble, France
³Flexcompute Inc, 130 Trapelo Road, Belmont, MA, 02478
⁴Dept. of Electrical & Computer Engineering, University of Wisconsin, Madison, WI, 53705
⁵Department of Applied Physics, Yale University, New Haven, Connecticut 06520
^⋆ yamilov@mst.edu
^† zongfu@flexcompute.com
^∗ hui.cao@yale.edu

1 Methods

1.1 Numerical method

Over the years, various techniques have been developed for large-scale 3D computations. One efficient method for simulating light scattering in large inhomogeneous media is based on Born series \citesupposnabrugge2016convergent, kruger2017solution. This method, however, is limited to small contrast in refractive index between the scatterers and the host medium. To achieve 3D localization, a large refractive-index contrast is needed to enhance light scattering, which makes it necessary to include many terms of the Born series. This greatly increases the computational complexity, worsening the convergence of such an iterative algorithm.

Another family of approaches for simulating wave transport in 3D aggregates of spherical particles is based on iterative multi-sphere scattering expansions, e.g., the generalized multi-particle-Mie (GMM) formalism \citesupp2001_Xu_GMM, 2007_Pellegrini_GMM, 2022_Egel_GMM, the multi-sphere T-matrix (MSTM) \citesupp2011_Mackowski_MSTM, and the fast multipole method (FMM) \citesuppgumerov2005computation. The latest, and the most advanced implementation, deployed on GPUs for speedup, is capable of simulating up to $10^{5}$ dielectric spheres \citesupp2017_Wiersma_3D_simulations, (?). But the state-of-the-art algorithm (CELES) has not been used for index contrast higher than $2$ . This is because the Mie resonances of high-index spheres can have very long lifetimes and strong internal fields, leading to an extremely slow convergence of the iterative algorithm or even its failure. Moreover, the spheres cannot touch or overlap, which would invalidate the Mie solution for isolated spheres. If two spheres are very close to each other, the near-field effects will excite high-order multipoles, causing slow or failed convergence. Therefore, these methods require a minimum distance between densely-packed spheres, which introduces spatial correlations of the scattering structures (?) that are known to affect Anderson localization \citesupp2008_Conti_AL_laser, (?).

The finite-difference time-domain (FDTD) method directly solves the Maxwell’s equations in space and time \citesupp2005_Taflove. It does not make any physical approximations, and is capable of simulating large-refractive-index, spatially-overlapping or touching particles with high filling fraction. The effects of spatial correlations of scattering structures that could contribute to Anderson localization \citesupp2008_Conti_AL_laser, (?), are avoided by completely-random placement of individual scatterers in our study, as shown in Sec. 1.4 below.

The ability of simulating refractive-index contrast of 3.5–10 in our study is the key to demonstrate that increasing index contrast will not precipitate Anderson localization. However, simulating such high refractive-index particles requires very fine spatial and temporal resolution, leading to an extremely long computational time. Previous FDTD simulations were limited to small 3D structures containing $10^{3}$ – $10^{4}$ particles (?), (?). Our hardware-accelerated implementation of the FDTD method (?) reduces the computational time by several orders of magnitude, allowing us to simulate a 3D system with $6\times 10^{6}$ scatterers in about 40 min.

The orders-of-magnitude reduction in computational time is essential to search for Anderson localization in 3D PEC and real-metal composites, because it allows to repeat the simulations many times with varying parameters like particle size, volume fraction, dielectric function and conductivity, lateral dimension and thickness of simulated systems, etc. Moreover, our time-domain simulation can extract the fields at more than $1000$ discrete frequencies from a single run, by Fourier transform. In contrast, the frequency-domain methods based on Born series and GMM/MSTM simulate only one frequency per run. Using such methods to simulate time-dependent transport, particularly at long delay time, requires repeating the calculations at many closely-spaced frequencies so that the Fourier transform will give the long-time behavior.

1.2 System geometry and dimension

Unless otherwise specified, our simulations are carried out on disordered systems in the slab geometry $L_{x}=L_{y}\gg~{}L$ (Fig. 1a). A plane wave with the electric field polarized along $x$ -axis is incident on the front surface of the slab. The bandwidth of a Gaussian pulse, $\Delta\lambda=$ 90 nm, is chosen such that $\ell_{t}$ remains nearly invariant with wavelength. Periodic boundary conditions are applied along $x$ and $y$ axes. To ensure that the periodicity would not affect the results presented in this manuscript, we make the transverse dimensions $L_{x},L_{y}$ of the simulated slabs much larger than the thickness $L$ . In the transmission simulations such as those in Fig. 1, the ratio is $L_{x}/L=L_{y}/L=3$ , so that any effect induced by periodic boundary conditions occurs on a length scale much larger than any relevant transport or localization scales. For the simulations of the transverse spreading (Figs. 4,S2,S12), we simulate much wider slabs having $L_{x}/L=L_{y}/L=10$ to avoid any edge effect.

The slab is sandwiched between homogeneous layers of refractive index $n_{B}$ and thickness $\Delta_{0}$ , which in turn are surrounded by perfectly matched layers of thickness $\Delta_{PML}$ . These thicknesses are chosen as $\Delta_{B}=\Delta_{PML}=2\lambda_{0},\lambda_{0},2\lambda_{0}$ for the slabs of dielectric spheres with refractive index $n=3.5,10$ and PEC spheres, respectively. For the dielectric slabs, $n_{0}$ is equal to the effective index $n_{\text{eff}}$ found by averaging the dielectric constant, while for PEC $n_{0}=1$ . $\lambda_{0}=650$ nm is the central wavelength of optical pulses in all simulations. The spatial discretization step of the FDTD algorithm is $\lambda_{0}/20$ for $n=3.5$ and PEC, and $\lambda_{0}/60$ for $n=10$ . This corresponded to time steps of $dt\simeq 56\times 10^{-18}$ s and $dt\simeq 19\times 10^{-18}$ s, respectively.

1.3 Numerical accuracy and scaling

Tidy3D solver (?) is an implementation of the standard FDTD method, which does not make any physical approximation or impose any constraint of the scattering structures. The remaining numerical issue in FDTD simulation is sufficient discretization of space and time \citesupp2005_Taflove. We have carefully tested spatio-temporal resolution of FDTD simulations to ensure consistency of our numerical results so that the conclusions of our study are robust and independent of the discretization. We present two examples of such tests below.

(i) One key evidence for the absence of AL in 3D dielectric random media is the exponential decay of transmitted flux $T(t)$ , shown in Fig. 1e of the main text. This result is obtained with the spatial grid size of $\lambda_{0}/20$ . When the spatial grid is reduced to $\lambda_{0}/40$ and the temporal step size is also halved, the exponential decay of $T(t)$ persists over 12 orders of magnitude in Fig. S1, which confirms the diffusive transport.

(ii) A tell-tale sign of AL in PEC composites is the arrest of transverse spreading of the transmitted beam, when an incident pulse is focused to a slab, as shown in Fig. 4c,d of the main text. This result is obtained with numerical resolution of $\lambda_{0}/20$ . To test the robustness of this result, we vary the resolution from $\lambda_{0}/10$ to $\lambda_{0}/40$ . Fig. S2a shows the result with $\lambda_{0}/40$ resolution: the transverse diameter of transmitted beam $d(t)$ saturates to a constant value in time, consistent with the result of $\lambda_{0}/20$ resolution in Fig. 4d. The asymptotic beam diameter $d_{\infty}$ is plotted versus the numerical resolution in Fig. S2b. It confirms the consistency of our numerical results: the arrest of transverse spreading persists in all simulations performed on the progressively finer meshes. Notably, the arrest is already seen in simulation with $\lambda_{0}/10$ resolution, despite a slight difference in the value of $d_{\infty}$ .

Since the Tidy3D implements the standard FDTD algorithm, the scaling of computing time with system size is the same as the standard FDTD method \citesupp2005_Taflove. More specifically, the computing time scales as $N_{x}\,N_{y}\,N_{z}\,N_{t}$ , where $N_{x}$ , $N_{y}$ , $N_{z}$ denote the number of spatial grid points in $x$ , $y$ , $z$ dimensions, $N_{t}$ is the total number of time steps. Typically the spatial grid size along $x$ , $y$ , $z$ is identical, and it is proportional to the time step.

We note that the finite discretization introduces slight anisotropy and dispersion to the simulated system. To mitigate such effects, the spatio-temporal resolution could be further increased according to the system size, which would affect the scaling of computing time \citesuppkruger2017solution. However, this is not necessary in our simulations of random ensembles of dielectric or metal spheres, because the statistical properties like scattering and transport mean free paths are barely modified. By varying the spatio-temporal grid size, we find the small effects caused by finite discretization are simply equivalent to a tiny change of the refractive index $n$ of the medium and/or a tiny change of the volume filling fraction $f$ . As shown in the manuscript, AL is absent in the dielectric systems for a broad range of $n$ ’s, whereas AL persists in a broad range of $f$ ’s in the PEC systems.

Our extensive tests of the numerical procedure and accuracy confirm the consistency and robustness of the results and conclusions in this manuscript.

1.4 Spatial correlation function

In order to avoid any residual spatial correlation of scattering structures, we adopt a uniform random distribution of scatterer centers. The structure factor obtained from Fourier transform of center positions is equal to unity \citesupp2018_Torquato_review. The spheres with center spacing smaller than their diameter will overlap in space, and the overlapping region is assigned the refractive index equal to that of an isolated sphere.

We calculate the spatial correlation function $C(\Delta)=\langle h(\boldsymbol{\rho})\,h(\boldsymbol{\rho}+{\bf\Delta})\rangle/\langle h(\boldsymbol{\rho})\rangle^{2}-1$ , where, $\boldsymbol{\rho}$ denotes the spatial position, $h(\boldsymbol{\rho})$ is a binary function equal to $1$ inside the dielectric/PEC/metal scatterers and $0$ outside, and $\langle...\rangle$ denotes averaging over $\boldsymbol{\rho}$ , direction of ${\bf\Delta}$ , and random ensemble. Fig. S3 shows the normalized $C(\Delta)/C(0)$ with $h(\boldsymbol{\rho})$ obtained from the discretized structures in the actual FDTD calculations of systems with particle radii $r=50$ nm and $100$ nm, and the PEC volume fractions of $f\sim 15\%$ and $\sim 50\%$ . In all cases, the spatial correlation extends over the range $\Delta$ of one particle diameter $2r$ , as expected for a random arrangement of spherical particles.

1.5 Scattering mean free path

Scattering mean free path $\ell_{s}$ measures the average distance traveled between two consecutive scattering events. Its value in Figs. 1c and 3a is extracted from the attenuation of the coherent component of the incident field. We simulate systems with dimension $L_{x}=100\lambda_{0},L_{y}=L=2\lambda_{0}$ for $n=3.5$ and PEC, $L_{x}=100\lambda_{0},L_{y}=L=\lambda_{0}$ in the $n=10$ case. The quantity $\langle E_{x}(x,y_{0},z;\lambda)\rangle_{x}$ is computed by averaging over the long dimension (along $x$ -axis) of the system for one particular cross section $y=y_{0}$ . Average co-polarized field amplitude $|\langle E_{x}(x,y_{0},z;\lambda)\rangle_{x}|$ decays exponentially at a rate $1/(2\ell_{s})$ , from which $\ell_{s}$ is obtained for each wavelength. The phase of $\langle E_{x}(x,y_{0},z;\lambda)\rangle_{x}$ gives $k_{\text{eff}}$ .

Figures S4, S5, S6 show the amplitude and phase of coherent field $\langle E_{x}(x,y_{0},z;\lambda)\rangle_{x}$ , as well as its real and imaginary parts, in random aggregates of dielectric spheres and PEC spheres. In the dielectric systems, spectral regions with strong attenuation, i.e. short scattering mean free path, are concentrated in the vicinity of Mie resonances of the constituent spherical particles. In contrast, in the PEC systems, scattering mean free path does not exhibit notable spectral features in Fig. 3a. In Figs. S6a,e, this can be seen from a constant attenuation rate of the field amplitude.

1.6 Transport mean free path

Transport mean free path $\ell_{t}$ corresponds to the average travel distance that is required to completely randomize the propagation direction. Transmittance of a continuous wave (CW) at wavelength $\lambda$ through a diffusive slab of thickness $L$ is $T(\lambda)=(5/3)\ell_{t}(\lambda)/[L+2z_{0}(\lambda)]$ , where $z_{0}(\lambda)$ is the extrapolation length (?). The value of $z_{0}(\lambda)$ can be estimated from the CW depth profile of internal intensity $I(z,\lambda)$ by linear extrapolation $I(z=L+z_{0},\lambda)=0$ (Fig. S9c). Then the value of $\ell_{t}(\lambda)$ is extracted from $T(\lambda)$ with known $z_{0}(\lambda)$ .

Typically, $z_{0}$ depends on the refractive-index mismatch between the slab and the surrounding medium (?). However, our choice of the refractive index $n_{\text{eff}}$ for the surrounding layers eliminates this mismatch for a dielectric slab, leading to $z_{0}=(2/3)\ell_{t}$ (?). Our numerical data suggest that this relation approximately holds for PEC composites embedded in air as well.

1.7 Diffusion coefficient

Diffusion coefficient $D=\ell_{t}v_{E}/3$ depends on both the transport mean free path $\ell_{t}$ and the energy velocity $v_{E}$ \citesupp1996_Lagendijk. The propagation of an optical pulse in a slab geometry makes it possible to directly extract $D$ from the decay of the transmittance, $T(t)\propto\exp[-t/t_{D}]$ , where $1/t_{D}=\pi^{2}D/(L+2z_{0})^{2}$ . In the localized slabs such as those in Fig. 2, the decay rate, $1/t_{D}$ , changes with time. To obtain $D(t)$ in Fig. 2b, we use an exponential fit within a time window $[t-2\tau_{D},t+2\tau_{D}]$ , where $\tau_{D}$ is the arrival time of the peak of transmitted flux.

1.8 Scaling functions

To obtain the scaling function $d\log(T)/d\log(L)$ for the PEC slabs, we compute the logarithmic derivative of CW transmittance $T$ with respect to slab thickness $L$ for a range of values of $k\ell_{s}$ . First, we use the wavelength dependence of $k\ell_{s}$ in Fig. 3a to map $k\ell_{s}$ to $\lambda$ for each $f$ . Next, we compute $T(\lambda)$ for two systems with different thicknesses $L_{1}=2\lambda_{0}$ and $L_{2}=5\lambda_{0}$ , where $\lambda_{0}=650$ nm denotes the center of the wavelength range of interest. Finally, we approximate $d\log(T)/d\log(L)$ by the finite difference $[\log(T_{2})-\log(T_{1})]/[\log(L_{2})-\log(L_{1})]$ for all $\lambda$ and all $f$ . After eliminating $\lambda$ , we obtain the dependence of $d\log(T)/d\log(L)$ on the Ioffe-Regel parameter $k\ell_{s}$ shown in Fig. 3c.

2 Absence of light localization in 3D dielectric systems with refractive-index contrast of $10$

This section reports the results of numerical simulation for a 3D slab of dielectric spheres with $n=10$ . For ease of comparison, we adjust the sphere radius to $r=32$ nm, so that the first Mie resonance occurs at a wavelength similar to that for $n=3.5$ sphere (compare Figs S4a and S5a). Scattering cross section of a single sphere of $n=10$ is exceedingly large (Fig. S7a), leading to strong dependent scattering already at small volume filling fractions, and limiting Ioffe-Regel parameter to $k_{\text{eff}}\,\ell_{s}\gtrsim 1$ in Fig. S7b. The minimal value of $k_{\text{eff}}\,\ell_{s}$ is almost identical for $f=2.5\%,\,5\%,\,9\%$ , and starts to increase for $f=18\%$ due to dependent scattering. The dependent scattering also results in less variation of the normalized transport mean free path $k_{\text{eff}}\,\ell_{t}$ with $\lambda$ in Fig. S7c, similar to the scattering mean free path in Fig. S7b.

Diffusive nature of wave propagation in $n=10$ dielectric systems can be seen from scaling of the CW transmittance with slab thickness $L$ . In Fig. S7d we plot transmittance multiplied by $L$ in the $f=2.5\%,\,5\%,\,9\%,\,18\%$ slabs with three different values of $L/\lambda_{0}=1,2,4$ . The largest slab thickness $L/\lambda_{0}=4$ corresponds to $L/\ell_{t}>36\gg 1$ for all $f$ . The overlap between the curves with different $L$ in the spectral range of the first Mie resonance is a direct manifestation of the diffusive scaling $T\propto 1/L$ .

Furthermore, we numerically calculate dynamic transmittance under pulsed excitation, $T(t)$ . For thick slabs ( $L/\ell_{t}\gg 1$ ), it does not exhibit deviation from the diffusive transport: at $t\gg\tau_{D}$ , the decay is still exponential (Fig. S7e).

The above results confirm light transport is diffusive in 3D dielectric random systems with $n=10$ spheres.

3 Difference between PEC and dielectric scattering

Besides the absence of longitudinal fields in PEC composites, additional differences with respect to dielectric media facilitate Anderson localization in PEC:

1.

Isolated PEC spheres have strong backscattering (?). It results in a highly anisotropic angular scattering pattern, making $\ell_{t}<\ell_{s}$ . This is very different from the dielectric spheres.
2.

Collective scattering resonances are created by two or more adjacent PEC spheres. The spatial arrangement of scatterers in this work lacks spatial correlations in sphere positions, creating a random distribution of distance between them. Even at the PEC volume fraction of $f=8\%$ , the gap in-between two PEC spheres can be much narrower than $\lambda=650$ nm, and the field intensity in the vicinity of the gap is strongly enhanced, as seen in Figs. S8a. Such an enhancement persists at $\lambda=1080$ nm in Figs. S8b, despite the wavelength is much greater than the particle radius $r$ = 50 nm, leading to a small size parameter $kr\simeq 0.3<1$ . Already for volume fraction of $f=8\%$ , a sufficiently broad distribution of gap sizes creates a large number of red-shifted resonances, enhancing the overall scattering at longer wavelengths in Fig. 3a, and making $\ell_{s}$ well below the prediction \citesupp1999_van_Rossum of the independent scattering approximation (ISA) $\ell_{s}=[\rho_{sca}(f)\,\sigma_{sca}(\lambda)]^{-1}$ , where $\rho_{sca}(f)$ is the number density of scatterers and $\sigma_{sca}(\lambda)$ is the scattering cross section of a single PEC sphere. With an increase of $f$ to $15\%$ in Figs. S8c,d, additional voids formed by three or more PEC spheres create even a larger variety of scattering resonances at different wavelengths rendering $\ell_{s}$ nearly constant with $\lambda$ .
3.

At high volume fraction $f$ of PEC scatterers, light propagation through air voids is suppressed. In contrast to a dielectric composite which becomes transparent when the dielectric volume fraction approaches 100%, a PEC composite with $f$ approaching 100% becomes a perfect mirror, as air voids no longer percolate through the system. In our simulations, Anderson localization takes place at PEC volume fractions well below the air percolation threshold (96.6%) \citesupp1981_Percolation_of_spheres,1984_Elam_Percolation_of_spheres. The air voids are still connected by narrow channels through which light may propagate. However, when the typical size of air voids between the PEC scatterers is on the order of the wavelength, light can no longer ‘squeeze’ through the narrow channels. Thus AL occurs at the wavelength comparable to or larger than the typical size $a$ of air voids between PEC scatterers.

For a given PEC volume fraction $f$ , air void size $a$ is proportional to PEC sphere radius $r$ . A random system with smaller spheres have a larger number of voids, but the voids themselves have smaller sizes, since the total volume of all voids is fixed by $f$ . The smaller voids make it harder for light to propagate through a system with smaller spheres, and localization will take place at lower $f$ . This prediction agrees to our numerical results of lower critical $f$ for AL with $r=50$ nm PEC spheres than $r=100$ nm spheres. The above argument is for fixed wavelength $\lambda$ . In a PEC system with given $r$ and $f$ , the typical void size $a$ is fixed, increasing $\lambda$ will make it more difficult for light to ‘squeeze’ through voids, facilitating AL at long wavelengths. Consequently, for a given $r$ , the localization condition of $a\sim\lambda$ is fulfilled at a smaller $f$ at longer $\lambda$ , which agrees to our numerical results.

4 Additional evidence for Anderson localization in PEC system

We have examined the intensity distributions inside 3D slabs of PEC spheres with $r=50$ nm in dynamical simulations. Fig. S9a shows the cross-section averaged depth profiles at very long delay times, $\langle I(x,y,z;t\rightarrow\infty)\rangle_{x,y}$ , for PEC filling fractions of $15\%$ and $48\%$ . In the diffusive slab of $f=15\%$ , the asymptotic depth profile matches the first eigenmode of the diffusion equation, $\sin[\pi(z+z_{0})/(L+2z_{0})]$ . In contrast, the cross-section averaged intensity inside the localized slab with $f=48\%$ exhibits much larger variation with depth $z$ . The faster decay towards the surfaces reflects stronger confinement of energy near the center of the slab.

Fluctuation of intensity or transmission coefficients is a powerful criterion of localization (?). Fig. S9b depicts spectral fluctuation of the cross-section average intensity $I(z,\lambda)=\langle I(x,y,z;\lambda)\rangle_{x,y}$ , within the wavelength range of 600–700 nm. The normalized variance
${\rm var}_{\lambda}[I(z,\lambda)]/\langle I(z,\lambda)\rangle_{\lambda}^{2}$ corresponds to the magnitude of long-range correlation $C_{2}$ \citesupp2014_Sarma_Cofz. In the PEC slab of $f=15\%$ , $C_{2}\ll 1$ , typical of a diffusive system. When $f$ increases to $48\%$ , $C_{2}\geq 1$ , as expected for localized systems\citesupp1999_van_Rossum,2004_Yamilov_passive_corr, (?).

Inside a diffusive system, the CW intensity decays linearly with the depth (?), and exhibits low fluctuations from realization to realization. This is confirmed in Fig. S9c by comparing $\langle I(z,\lambda)\rangle_{\lambda}$ and $\exp[\langle\log[I(z,\lambda)]\rangle_{\lambda}]$ , which indeed agree well in the PEC slab of $f=15\%$ . In contrast, in the PEC slab of $f=48\%$ , these two quantities are markedly different owning to the fact that strong intensity fluctuations lead to log-normal distribution for localized systems \citesupp2000_Mirlin. The depth profile exhibits a roughly exponential decay, as a result of AL.

To further confirm AL in PEC slabs, we repeat the calculations reported in Fig. 3c,d for another system with larger PEC spheres ( $r=100$ nm). Figure S10 demonstrates that transmittance and conductance exhibit similar scaling behaviors as for the PEC slabs of $r$ = 50 nm in Fig. 3c,d. Although the diffusion-localization transition occurs at a different volume filling fraction $f$ of PEC spheres, the scaling behaviors with respect to the Ioffe-Regel parameter $k\,\ell_{s}$ and to the dimensionless conductance $g$ are reproduced.

5 Anderson localization in real metals

To demonstrate the possibility of AL in realistic metallic systems, we simulate real metals using parameters reported in literature. The dielectric constant of a metal is related to the conductivity $\sigma(\omega)$ via $\epsilon(\omega)=1+i\sigma(\omega)/\omega\epsilon_{0}$ , where $\epsilon_{0}$ is electric permittivity of vacuum. In the Drude model, the conductivity of a metal is given by $\sigma(\omega)=\sigma_{0}/(1-i\omega\tau)$ , where $\sigma_{0}$ is the static conductivity and $\tau$ is the relaxation time.

5.1 Microwave regime

At microwave frequencies, $\omega\ll 1/\tau$ and $\epsilon(\omega)\simeq 4\pi\sigma_{0}i/\omega$ . The penetration depth of electric field into the metal is characterized by the skin depth $\sim(\lambda_{0}/2\pi)\,(\omega/2\pi\sigma_{0})^{1/2}$ . Common metals like silver, aluminum, and copper have low loss and large conductivity at microwave frequencies, and their skin depth is much shorter than the wavelength. For example, at $20$ GHz where the wavelength is $1.5$ cm, the skin depth is less than $1\ \mu$ m \citesupp2017_Yi_Plasmonics_book.

To simulate microwave transport in 3D aggregates of metallic spheres, we generate random arrangements of overlapping spheres using the same procedure as for the PEC systems. The sphere diameter of $0.56$ cm is four orders of magnitude larger than the skin depth of crystalline metals like silver, aluminum, copper. Thus the microwave barely penetrates into the metallic scatterers, and the absorption loss is negligible. When we use the conductivity $\sigma_{0}=3.8\times 10^{7}$ ${\rm\Omega^{-1}/m}$ of crystalline aluminum reported in literature \citesupp2017_Yi_Plasmonics_book, the simulation results barely deviate from those for PEC.

However, experimentally fabricated metallic structures have additional loss due to polycrystallinity, surface defects, oxide layers, etc. To take these into account, we lower the conductivity $\sigma_{0}$ so that the simulated transport properties match the experimental values in Ref. (?). Specifically, we simulate dynamic transmittance in a diffusive slab of aluminum spheres with same volume fraction $f=35\%$ and slab thickness $L$ = 6 cm as in the experiment. When the conductivity is reduced to $\sigma_{0}=3.8\times 10^{4}$ ${\rm\Omega^{-1}/m}$ , our simulated diffusion coefficient $D=1.9\times 10^{9}$ cm²/s and absorption coefficient $\alpha=(D\tau_{a})^{-1}\simeq 0.2$ cm^-1 both agree to the experimental values at 20 GHz. Using these realistic parameters, we simulate 3D metallic composites with high volume fraction of $f=60\%$ . As shown in Fig. S11a, the apparent diffusion coefficient, obtained from the temporal decay of the transmittance without taking absorption into account, decreases in time as $1/t$ , similarly to its behavior in PEC, until it reaches the value set by the absorption. The simulated field intensity distribution inside the system, in Fig. S11b, provides evidence of spatial confinement of microwave, similar to the result for PEC in Fig. 2e. The most conclusive evidence of AL is the arrest of transverse spreading of transmitted wave field. It is insensitive to absorption and is shown in the main text in Fig. 4e.

We believe the tell-tale experimental sign of AL is the arrest of transverse spreading of transmitted wave with a focused incident pulse. The potential pitfalls are background signals or possible emission by the metal upon microwave excitation. However, these signals are typically broadband and incoherent with the incident wave. If the incident microwave has a narrow spectral band, the transmitted field may be measured coherently, e.g., using the homodyne technique common for microwaves, and those background and incoherent signals will be discarded. Hence, we propose an experiment with a frequency-tunable narrowband microwave source: focus the incident microwave to the front surface of a slab and scan the frequency, perform an interferometric measurement of transmitted field distribution near the slab back surface at each frequency (using a local oscillator), finally Fourier transform the spectral fields to reconstruct the temporal evolution of transmitted field for an incident short pulse in order to observe the arrest of transverse spreading.

5.2 Optical regime

Here we investigate the possibility of AL of visible light in 3D metallic nanostructures. Even for low-loss metals like gold and silver, the skin depth is $\sim 25$ nm at $\lambda=650$ nm \citesupp2017_Yi_Plasmonics_book, comparable to the nanoparticle diameter of 100–200 nm. A significant penetration of light field inside the metal particles makes them deviate from the PEC, which expels the fields completely. Another consequence of the penetration is a notable loss. We simulate silver, which is widely used in nano-plasmonics and meta-materials, using realistic parameters reported in the literature. We adopt the Drude model of $\epsilon(\omega)$ with $\sigma_{0}\simeq 6.1\times 10^{7}$ ${\rm\Omega^{-1}/m}$ and $\tau\simeq 3.7\times 10^{-14}$ s from Ref. \citesupp1985_Ordal_Drude_model_parameters. Despite absorption and deviation from PEC, we still observe the arrest of transverse spreading in Fig. S12.

Compared to 2D nanostructures, 3D nanoporous metals have a much larger (internal) surface area, leading to a wide range of applications in photo-catalysis \citesupplinic2011plasmonic, optical sensing \citesuppchen2009nanoporous, energy conversion and storage \citesuppmascaretti2019plasmon. Metallic nanostructures have been widely explored to enhance Raman scattering, second-harmonic generation, etc., because they can produce ‘hot spots’ (giant local fields) to boost optical nonlinearities. Also metallic nanoparticles have been used for random lasing, as their strong scattering of light improves optical feedback \citesuppwang2017nanolasers, gomes2021recent. Our simulation results suggest the possibility of AL in 3D metallic nanostructures at optical frequencies. Light localization in such structures will have a significant impact on optical nonlinear, lasing, photochemical processes and related applications.

\bibliographystylesupp

Nature \bibliographysuppsupp.bbl