# Anderson localization of electromagnetic waves in three dimensions

Alexey Yamilov, Sergey Skipetrov, Tyler Hughes, Momchil Minkov, Zongfu Yu, Hui Cao

► To cite this version:

Alexey Yamilov, Sergey Skipetrov, Tyler Hughes, Momchil Minkov, Zongfu Yu, et al., Anderson localization of electromagnetic waves in three dimensions. Nature Physics, 2023, 19 (9), pp.1308-1313.
10.1038/s41567-023-02091-7. hal-04134647

HAL Id: hal-04134647

https://hal.science/hal-04134647

Submitted on 20 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Anderson localization of electromagnetic waves in
three dimensions

Alexey Yamilov¹*, Sergey E. Skipetrov², Tyler W. Hughes³, Momchil Minkov³,
Zongfu Yu³,⁴,†, Hui Cao⁵*

¹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

September 16, 2022

Anderson localization marks a halt of diffusive wave propagation in disor-
dered 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 numer-
ical simulations of light transport in fully disordered 3D systems with un-
precedented 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 dielec-
tric spheres with a refractive index contrast up to 10. Our work opens a wide
computing hardware has brought orders of magnitude speed-up of numerical calculation (36, 37). 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 (38), the interplay between disorder and topological phenomena may be explored beyond the limit of weak disorder in low-dimensional systems (39). Also in cavity quantum electrodynamics with Anderson-localized modes (40), going 3D will avoid the out-of-plane loss inherent for 2D systems and cover the full angular range of propagation directions (41). In addition to fundamental studies, disorder and scattering has been harnessed for various photonic device applications, but mostly with diffuse waves (42). 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 (43). Rate of this exponential decay is equal to $1/\tau_D$, which is directly related (43) 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 (8, 23).

Previous studies (26, 27) 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 E(\mathbf{r}) = 0$ does not follow from the Gauss's law $\nabla \cdot [\epsilon(\mathbf{r})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) (44) 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 (8).

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 µm × 10 µm × 3.3 µ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 (45). The non-exponential decay of $T(t)$ and the time-dependence of $D$ are the signatures of AL (13, 45). 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 ultimately related to their spatial confinement (3, 9, 46). 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 (47).

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 $kl_s$. We compute it using a procedure similar to that in dielectrics (SI). Even at the volume fraction of $f = 8\%$, $l_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, $l_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 \approx 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 (48) 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$ (49) 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})$ (49). In diffusive regime $g > 1$, $\beta(g) \to 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 (13, 33, 50). A pulse centered at $\lambda = 650$ nm with a bandwidth of 90 nm is focused to a small spot of area $\approx 0.5$ $\mu$m$^2$ at the front surface of a wide 3D slab of dimensions 33 $\mu$m×33 $\mu$m×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[PR(t)/\pi]^{1/2}$, where
$$PR(t) = [\iint I(x, y, L; t)dxdy]^2 / \iint 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_{eff} l_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 (29). 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 (22, 23, 31–33). 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 (29). In the SI, we propose a realistic microwave experiment that avoids experimental pitfalls (9) 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 *Ultra-cold 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).
<figure><img src="image_2.png" /><figcaption>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^2$ 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 $\mu m^2/s$ at $f \approx 30\%$.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: **Anderson localization of light in 3D disordered PEC.** a, Transmittance $T(t)$ of an optical pulse through a 3D slab (10 µm × 10 µm × 3.3 µm) of randomly packed PEC spheres with radius $r = 50$ nm and volume filling fractions $f$ from 8% to 48%. b, Time-resolved diffusion coefficient $D(t)$ extracted from the decay rate of $T(t)$ in panel a, decreasing with time as $1/t$ at high $f$. c, Short-time $D$ (dots) and the interval of variation of $D$ with time (bars) at different $f$. d,e, Continuous-wave transmittance spectrum $T(\lambda)$ in diffusive (d, $f = 15\%$, blue line) and localized (e, $f = 48\%$, red line) PEC slabs. Color map: depth profile of average intensity $\langle I(x, y_0, z; \lambda) \rangle_x$ inside the slab at different wavelengths, highlighting the localized and necklace-like states for $f = 48\%$.</figcaption></figure>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

Figure 3: **Transition from diffusion to Anderson localization in 3D disordered PEC.** a, Scattering mean free path $l_s$ for PEC volume filling fractions $f$ from 8% to 48%. $l_s$ is nearly flat over broad spectral range. b, Spectral dependence of Ioffe-Regel parameter $kl_s$, exhibiting $1/\lambda$ dependence (dashed lines). c, Scaling of the continuous-wave (CW) transmittance $T$ with slab thickness $L$ versus Ioffe-Regel parameter $kl_s$, revealing diffusion-localization transition at $kl_s \sim 1$. Blue dashed line denotes diffusive scaling $T \propto L^{-1}$; red dashed line marks $kl_s = 1$. d, Single-parameter scaling of dimensionless conductance $g$ for diffusion-localization transition (solid black line), in agreement with numerical data for 6 PEC filling fractions. Blue and red dashed lines denote diffusive and localized scalings $g \propto L$ and $g \propto \exp(-L/\xi)$ respectively, where $\xi$ is the localization length.
<figure><img src="image_8.png" /><figcaption>Figure 4: Arrest of transverse spreading of transmitted beam in 3D localized PEC systems.
a, Schematic of transverse spreading of a tightly-focused pulse propagating through a diffusive slab of cross section 33 µm × 33 µm and thickness L = 3.3 µm. b, 2D intensity distribution at the output surface (normalized to the maximum) for different delay times, showing a lateral expansion of the beam in the diffusive slab of PEC filling fraction f = 15%. c, Absence of transverse spreading for f = 48%, due to Anderson localization. d, Lateral diameter of the transmitted beam d(t) increases as $\sqrt{t}$ in the diffusive slab (blue dots), but it saturates to a constant value in the localized slab (red dots). e, Same as d but for a slab (L = 6 cm) of aluminum spheres (r = 0.28 cm, f = 35% and 60%) with realistic conductivity $\sigma_0 = 3.8 \times 10^4$ Ω⁻¹/m at microwave frequency ~ 20 GHz.</figcaption></figure>
## 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<fnref n="1" />*, Sergey E. Skipetrov<fnref n="2" />, Tyler W. Hughes<fnref n="3" />, Momchil Minkov<fnref n="3" />,

Zongfu Yu<fnref n="3" />,<fnref n="4" />,†, Hui Cao<fnref n="5" />*

# 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 (S1, S2). 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-

<footnote n="1">¹Physics Department, Missouri University of Science & Technology, Rolla, Missouri 65409</footnote>

<footnote n="2">²Univ. Grenoble Alpes, CNRS, LPMMC, 38000 Grenoble, France</footnote>

<footnote n="3">³Flexcompute Inc, 130 Trapelo Road, Belmont, MA, 02478</footnote>

<footnote n="4">⁴Dept. of Electrical & Computer Engineering, University of Wisconsin, Madison, WI, 53705</footnote>

<footnote n="5">⁵Department of Applied Physics, Yale University, New Haven, Connecticut 06520</footnote>

<footnote>* yamilov@mst.edu</footnote>

<footnote>† zongfu@flexcompute.com</footnote>

<footnote>* hui.cao@yale.edu</footnote>
particle-Mie (GMM) formalism (S3, S4, S5), the multi-sphere T-matrix (MSTM) (S6), and the fast multipole method (FMM) (S7). The latest, and the most advanced implementation, deployed on GPUs for speedup, is capable of simulating up to $10^5$ dielectric spheres (S8), (35). 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 (35) that are known to affect Anderson localization (S9), (21).

The finite-difference time-domain (FDTD) method directly solves the Maxwell's equations in space and time (S10). 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 (S9), (21), 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 (34), (21). Our hardware-accelerated implementation of the FDTD method (36) 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 Ander-
son localization in 3D PEC and real-metal composites, because it allows to repeat the simula-
tions 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 thick-
nesses 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 (36) 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 (S10). 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 (S10). 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 (S2). 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 (S11). 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(\rho) h(\rho + \Delta) \rangle / \langle h(\rho) \rangle^2 - 1$, where, $\rho$ denotes the spatial position, $h(\rho)$ is a binary function equal to 1 inside the dielectric/PEC/metal scatterers and 0 outside, and $\langle \dots \rangle$ denotes averaging over $\rho$, direction of $\Delta$, and random ensemble. Fig. S3 shows the normalized $C(\Delta)/C(0)$ with $h(\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 $l_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/(2l_s)$, from which $l_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 (43). 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 (43). 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$ (43). 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$ (S12). 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 (44). 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 \sim 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 (S13) 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 $l_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%) (S14, S15). 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 \to \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 (9). 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 $\text{var}_\lambda[I(z, \lambda)]/\langle I(z, \lambda) \rangle_\lambda^2$ corresponds to the magnitude of long-range correlation $C_2(S16)$. In the PEC slab of $f = 15\%$, $C_2 \ll 1$, typical of a diffusive system. When $f$ increases to 48%, $C_2 \ge 1$, as expected for localized systems (S13, S17), (9).

Inside a diffusive system, the CW intensity decays linearly with the depth (43), 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 owing to the fact that strong intensity fluctuations lead to log-normal distribution for localized systems (S18). 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* *l*<sub>s</sub>
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
σ(ω) via ϵ(ω) = 1+iσ(ω)/ωϵ₀, where ϵ₀ is electric permittivity of vacuum. In the Drude model,
the conductivity of a metal is given by σ(ω) = σ₀/(1-iωτ), where σ₀ is the static conductivity
and τ is the relaxation time.

**5.1 Microwave regime**

At microwave frequencies, $\omega \ll 1/\tau$ and $\epsilon(\omega) \approx 4\pi\sigma_0i/\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 (S19).

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 σ₀ = 3.8 × 10⁷
Ω⁻¹/m of crystalline aluminum reported in literature (S19), the simulation results barely deviate
from those for PEC.

However, experimentally fabricated metallic structures have additional loss due to polycrys-
tallinity, 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. (29). 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 \Omega^{-1}/m$, our simulated diffusion coefficient $D = 1.9 \times 10^9$ cm$^2$/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 ~ 25 nm at $λ = 650$ nm (S19), 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 \Omega^{-1}/m$ and $\tau \simeq 3.7 \times 10^{-14}$ s from Ref. (S20). 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 (S21), optical sensing (S22), energy conversion and storage (S23). 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 (S24, S25). 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.
**Supplementary References**

[S1] Osnabrugge, G., Leedumrongwatthanakun, S. & Vellekoop, I. M. A convergent born series for solving the inhomogeneous helmholtz equation in arbitrarily large media. *Journal of computational physics* **322**, 113–124 (2016).

[S2] Krüger, B., Brenner, T. & Kienle, A. Solution of the inhomogeneous maxwell's equations using a born series. *Optics express* **25**, 25165–25182 (2017).

[S3] Xu, Y. L. & Gustafson, B. A generalized multiparticle mie-solution: further experimental verification. *J. Quant. Spectr. Rad. Transf.* **70**, 395–419 (2001). Light Scattering by Non-Spherical Particles.

[S4] Pellegrini, G., Mattei, G., Bello, V. & Mazzoldi, P. Interacting metal nanoparticles: Optical properties from nanoparticle dimers to core-satellite systems. *Mat. Sci. Eng. C* **27**, 1347–1350 (2007).

[S5] Egel, A. *et al*. Smuthi: A python package for the simulation of light scattering by multiple particles near or between planar interfaces. *J. Quant. Spectr. Rad. Transf.* **273**, 107846 (2021).

[S6] Mackowski, D. & Mishchenko, M. A multiple sphere t-matrix fortran code for use on parallel computer clusters. *J. Quant. Spectr. Rad. Transf.* **112**, 2182–2192 (2011).

[S7] Gumerov, N. A. & Duraiswami, R. Computation of scattering from clusters of spheres using the fast multipole method. *The Journal of the Acoustical Society of America* **117**, 1744–1761 (2005).
[S8] Egel, A., Pattelli, L., Mazzamuto, G., Wiersma, D. S. & Lemmer, U. Celes: Cuda-accelerated simulation of electromagnetic scattering by large ensembles of spheres. *J. Quantitative Spectroscopy and Radiative Transfer* **199**, 103–110 (2017).

[S9] Conti, C. & Fratalocchi, A. Dynamic light diffusion, three-dimensional Anderson localization and lasing in inverted opals. *Nat. Phys.* **4**, 794–798 (2008).

[S10] Taflove, A., Hagness, S. C. & Piket-May, M. *Computational electromagnetics: the finite-difference time-domain method*, 629–670 (Elsevier, 2005).

[S11] Torquato, S. Hyperuniform states of matter. *Phys. Rep.* **745**, 1–95 (2018). Hyperuniform States of Matter.

[S12] Lagendijk, A. & van Tiggelen, B. A. Resonant multiple scattering of light. *Phys. Rep.* **270**, 143–215 (1996).

[S13] van Rossum, M. C. & Nieuwenhuizen, T. M. Multiple scattering of classical waves: microscopy, mesoscopy, and diffusion. *Rev. Mod. Phys.* **71**, 313–371 (1999).

[S14] Kertész, J. Percolation of holes between overlapping spheres : Monte carlo calculation of the critical volume fraction. *J. Physique Lett.* **42**, 393–395 (1981).

[S15] Elam, W. T., Kerstein, A. R. & Rehr, J. J. Critical properties of the void percolation problem for spheres. *Phys. Rev. Lett.* **52**, 1516–1519 (1984).

[S16] Sarma, R., Yamilov, A., Neupane, P., Shapiro, B. & Cao, H. Probing long-range intensity correlations inside disordered photonic nanostructures. *Phys. Rev. B* **90**, 014203 (2014).

[S17] Chang, S. H., Taflove, A., Yamilov, A., Burin, A. & Cao, H. Numerical study of light correlations in a random medium close to the Anderson localization threshold. *Opt. Lett.* **29**, 917–919 (2004).
[S18] Mirlin, A. Statistics of energy levels and eigen-functions in disordered systems. *Phys. Rep.* **326**, 259–382 (2000).

[S19] Li, Y. (ed.) *Plasmonic Optics: Theory and Applications* (SPIE press, 2017).

[S20] Ordal, M. A., Bell, R. J., Alexander, R. W., Long, L. L. & Querry, M. R. Optical properties of fourteen metals in the infrared and far infrared: Al, co, cu, au, fe, pb, mo, ni, pd, pt, ag, ti, v, and w. *Appl. Opt.* **24**, 4493–4499 (1985).

[S21] Linic, S., Christopher, P. & Ingram, D. B. Plasmonic-metal nanostructures for efficient conversion of solar to chemical energy. *Nature materials* **10**, 911–921 (2011).

[S22] Chen, L.-Y., Yu, J.-S., Fujita, T. & Chen, M.-W. Nanoporous copper with tunable nanoporosity for sers applications. *Advanced Functional Materials* **19**, 1221–1226 (2009).

[S23] Mascaretti, L. et al. Plasmon-enhanced photoelectrochemical water splitting for efficient renewable energy storage. *Advanced Materials* **31**, 1805513 (2019).

[S24] Wang, Z., Meng, X., Kildishev, A. V., Boltasseva, A. & Shalaev, V. M. Nanolasers enabled by metallic nanoparticles: from spasers to random lasers. *Laser & Photonics Reviews* **11**, 1700212 (2017).

[S25] Gomes, A. S., Moura, A. L., de Araújo, C. B. & Raposo, E. P. Recent advances and applications of random lasers and random fiber lasers. *Progress in Quantum Electronics* **78**, 100343 (2021).
<figure><img src="image_9.png" /><figcaption>Figure S1: Confirmation of the absence of Anderson localization in random dielectric media with refractive-index contrast of 3.5 with increased numerical resolution. Transmittance $T(t)$ of an optical pulse through a 3D slab of thickness $L = 3.3$ µm, filled with dielectric spheres at random uncorrelated positions (radius $r = 100$ nm, refractive index $n = 3.5$, volume filling fraction $f = 38\%$) in air. Numerical simulation is performed with the spatio-temporal discretization corresponding to $\lambda_0/40$, which is twice finer than that used to produce Fig. 1e in the main text. Pure exponential decay of $T(t)$ in time, over 12 orders of magnitude, is a hallmark of diffusive transport in the dielectric random system.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure S2: Robustness of the arrest of transverse spreading of transmitted beam in localized PEC systems against numerical discretization. The simulation is schematically depicted in Fig. 4a of the main text, and the structure parameters (sphere radius $r = 50$ nm, volume filling fraction $f = 48\%$) are identical to those in Fig. 4d. a, Transverse diameter of the transmitted beam $d(t)$, obtained from numerical simulation (slab thickness $L = 1.3~\mu$m) with $\lambda_0/40$ resolution, exhibits a saturation consistent with the result with $\lambda_0/20$ resolution in Fig. 4d. b, Arrest of the transverse spreading is observed in all simulations performed with numerical resolution from $\lambda_0/10$ to $\lambda_0/40$ with a consistent asymptotic transverse beam diameter $d_\infty$. Symbols and error-bars represent the mean value and standard deviation of $d(t)$ in the time interval 1 ps $< t < 2$ ps (gray area in panel a).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure S3: **Spatial correlation in PEC composites.** Normalized spatial correlation function $C(\Delta)/C(0)$ for random aggregates of overlapping spheres with radii $r = 50$ nm, 100 nm and volume filling fractions $f \simeq 15\%, 50\%$. Spatial correlation vanishes beyond one sphere diameter $2r$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure S4: Extinction of coherent field in a dielectric slab of n = 3.5 spheres. a, Scattering cross section of a single dielectric sphere of radius r = 100 nm and refractive index n = 3.5 in air. The spectral peaks correspond to Mie resonances. b,c,d,e, amplitude (b), phase (c), real (d) and imaginary (e) parts of coherent field $\langle E_x(x, y_0, z; \lambda) \rangle_x$ versus depth z and wavelength λ for dielectric filling fraction f = 29%. White lines represent z = mλ$_{eff}$, where m is an integer and λ$_{eff}$ = λ$_0$/n$_{eff}$. The coherent field amplitude experiences an enhanced extinction in the wavelength range 500–700 nm due to hybridized Mie resonances.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure S5: **Extinction of coherent field in a dielectric slab of n = 10 spheres.** a, Scattering cross section of a single dielectric sphere of radius r = 32 nm and refractive index n = 10 in air. The peak represents the first Mie resonance. b,c,d,e, amplitude (b), phase (c), real (d) and imaginary (e) parts of coherent field $\langle E_x(x, y_0, z; \lambda) \rangle_x$ versus depth z and wavelength λ for dielectric filling fraction f = 29%. White lines represent z = mλ$_{eff}$, where m is an integer and λ$_{eff}$ = λ$_0$/n$_{eff}$. The coherent field amplitude experiences an enhanced extinction in the spectral range 630–660 nm in the vicinity of the first Mie resonance.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure S6: Comparison of coherent fields in diffusive and localized slabs of PEC spheres. Amplitude (a,e), phase (b,f), real (c,g) and imaginary (d,h) parts of coherent field $\langle E_x(x, y_0, z; \lambda) \rangle_x$ versus depth $z$ and wavelength $\lambda$ for PEC filling fraction $f = 15\%$ (a,b,c,d) and $48\%$ (e,f,g,h). White lines represent $z = m\lambda$, where $m$ is an integer. In the localized slab ($f = 48\%$), extinction of coherent field is much stronger and nearly independent of wavelength.</figcaption></figure>
<figure><img src="image_15.png" /></figure>

Figure S7: Resonant scattering and diffusion in random ensembles of n = 10 dielectric spheres. a, Scattering cross section of a single dielectric sphere of radius r = 32 nm and refractive index n = 10. The vertical dotted lines mark spectral width of the excitation pulse in panel e. b, Resonantly enhanced scattering mean free path ℓs in the spectral vicinity of the first Mie resonance of single-particle scattering. The dielectric filling fraction f = 2.5% (blue), f = 5% (red), 9% (green) and 18% (brown). The minimum Ioffe-Regel parameter koffels ~ 2. The horizontal dashed line marks the Ioffe-Regel criterion koffels = 1 for 3D localization. c, Normalized transport mean free path koffels in the spectral vicinity of the first Mie resonance, reflecting saturation by dependent scattering between f = 2.5% and 18%. d, CW transmittance T(λ) for three slab thicknesses L/λ₀ = 1, 2, 4 (solid, dashed, and dotted lines). T(λ) multiplied by L/λ₀ remains nearly independent of L at wavelengths around minimum transmission, reflecting the diffusive scaling of T ∝ 1/L. Results are color-coded as in panels b,c. For f = 5%, 9% and 18%, the curves are shifted up vertically by one, two and three decades for legibility. e, Transmittance of 3D slabs with thickness L = 2 × λ₀ for pulsed excitation, showing an exponential decay (dashed lines) in time at long delay. The decay rate is determined by the minimum diffusion coefficient within the pulse bandwidth. Legend shows values of L/ℓt in each system. For f = 5%, 9% and 18%, the curves are shifted down vertically by two, four and six decades for legibility.
<figure><img src="image_16.png" /><figcaption>Figure S8: **Spatial intensity distribution inside disordered PEC system.** Field intensity at a $(y, z)$ cross-section inside a random ensemble of PEC spheres with radius $r = 50$ nm. PEC filling fraction is $f = 8\%$ in panels a,b and $f = 15\%$ in panels c,d. $\lambda = 650$ nm in panels a,c and $1080$ nm in panels b,d. Bottom row shows a magnified view of each distribution in the top row, revealing strong intensity enhancement in air voids between PEC particles due to coupled resonances. The size of the voids can be much smaller than wavelength $\lambda$ indicated by the white scale bar.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

Figure S9: **Additional evidence for Anderson localization in disordered PEC.** a, Depth profile $\langle I(x, y, z; t \to \infty) \rangle_{x,y}$ inside the PEC slab of $f = 15\%$ (blue solid line) and $48\%$ (red solid line). The lowest-order diffusive mode profile (black dashed line) matches that of $f = 15\%$ (blue). b, Normalized variance of cross-section integrated intensity, showing two orders of magnitude enhancement of intensity fluctuations for $f = 48\%$ over $f = 15\%$. c,d, CW intensity depth profiles $\langle I(z, \lambda) \rangle_\lambda$ (solid lines) and $\exp[\langle\log[I(z, \lambda)]\rangle_\lambda]$ (dashed lines), compared to linear decay for $f = 15\%$ and exponential decay for $f = 48\%$ (black dashed lines).
<figure><img src="image_18.png" /><figcaption>Figure S10: **Diffusion-localization transition in another PEC system with larger spheres.**
The scaling behavior reported in Fig. 3b,c is reproduced for a system of larger PEC spheres with radius $r = 100$ nm, despite the diffusion-localization transition occurs at a different volume filling fraction $f$ of PEC spheres. In a, blue dashed line denotes diffusive scaling $T \propto L^{-1}$; red dashed line marks $k\ell_s = 1$. In b, blue and red dashed lines denote diffusive and localized scaling $g \propto L$ and $g \propto \exp(-L/\xi)$ respectively, where $\xi$ is the localization length.</figcaption></figure>
<figure><img src="image_19.png" /></figure>

**Figure S11: Localization of microwaves in a 3D real-metal composite.** A random aggregate of aluminum spheres with radius $r = 0.28$ cm, realistic conductivity $\sigma_0 = 3.8 \times 10^4$ $\Omega^{-1}/m$, and volume filling fraction $f = 60\%$ localizes microwave of frequency around 20 GHz in a slab of thickness $L = 6$ cm. a, Time-resolved apparent diffusion coefficient $D(t)$, extracted from the temporal decay rate of $T(t)$, decreases as $1/t$, due to Anderson localization. Eventually it saturates to the value set by metal absorption. b, Wavelength-resolved transmittance $T(\lambda)$ (red line) exhibits fluctuations. Color map: depth profile of average intensity $\langle I(x, y_0, z; \lambda) \rangle_x$ inside the slab at different wavelengths, highlighting spatially localized and necklace-like states.
<figure><img src="image_20.png" /><figcaption>Figure S12: Arrest of transverse spreading of visible light in 3D random aggregates of silver particles. Transverse diameter of the transmitted light $d(t)$, with a tightly-focused incident pulse of center wavelength $\lambda_0 = 650$ nm, increases as $\sqrt{t}$ in the diffusive slab ($f = 15\%$, blue dots), but it saturates to a constant value in the localized slab ($f = 48\%$, red dots). The slab thickness $L = 2$ µm. Silver spheres have radius $r = 50$ nm, $\sigma_0 \simeq 6.1 \times 10^7$ Ω$^{-1}$/m and $\tau \simeq 3.7 \times 10^{-14}$ s.</figcaption></figure>