## Abstract

Deep-subwavelength all-dielectric composite materials are believed to tightly obey the Maxwell Garnett effective medium theory. Here, we demonstrate that the Maxwell Garnett theory could break down due to evanescent fields in deep-subwavelength dielectric structures. By using two- and three-dimensional dielectric composite materials with inhomogeneities at a scale of $\lambda /100$, we show that local evanescent fields generally occur near the dielectric inhomogeneities. When tiny absorptive constituents are placed there, the absorption and transmission of the whole composite will show strong dependence on the positions of the absorptive constituents. The Maxwell Garnett theory fails to predict such position-dependent characteristics because it averages out the evanescent fields. By taking the distribution of the evanescent fields into consideration, we have made a correction to the Maxwell Garnett theory so that the position-dependent characteristics become predictable. We reveal not only the breakdown of the Maxwell Garnett theory, but also a unique phenomenon of “invisible” loss induced by the prohibition of electric fields at deep-subwavelength scales. We believe our work promises a route to control the macroscopic properties of composite materials without changing their composition, which is beyond the traditional Maxwell Garnett theory.

© 2021 Chinese Laser Press

## 1. INTRODUCTION

Controlling the optical properties of photonic materials is a subject of great importance. One important approach to construct advanced optical materials is to build composite materials by intermixing two or more homogeneous constituents at deep-subwavelength scales. Generally, such a composite material can be homogenized as a continuous effective medium with uniform properties [Fig. 1(a)]. In 1904, Maxwell Garnett developed a simple, but immensely successful, effective medium theory (EMT) [1]. Fields within the unit cell are considered local, and material properties are determined by the polarization and magnetization vectors. In this way, the Maxwell Garnett EMT gives the effective permittivity (or permeability) of the composite in terms of the permittivities (or permeabilities) and filling ratios of the individual constituents of the composite material. In this scenario, the macroscopic optical properties (e.g., reflection, transmission, and absorption) of the composite material are insensitive to the details of the constituents (e.g., the position of each inclusion), as they are averaged out in the effective medium description [2–5].

Local homogenization is insufficient for the correct description of wave behaviors in special composite structures involving extremely large wave vectors or surface wave resonances even in the deep-subwavelength scale [6–13]. In such cases, the effective parameters become nonlocal or spatially dispersive (i.e., dependent on wave vectors) [14,15]. For example, strong nonlocality can be induced by surface plasmons in metal–dielectric structures even at the deep-subwavelength scale, leading to unusual effects that include additional modes [8] and parabolic dispersions [7,10].

On the other hand, deep-subwavelength all-dielectric composite materials, where surface wave resonances are not supported, are generally believed to tightly obey the Maxwell Garnett EMT. Interestingly, Sheinfux *et al.* [16] very recently showed the breakdown of EMT in deep-subwavelength all-dielectric multilayers, which has been numerically and experimentally demonstrated [17–27]. They found that the transmission through the multilayer structure depends strongly on nanoscale variations at the vicinity of the effective medium’s critical angle for total internal reflection. In this circumstance, the transmission spectra of the actual multilayer and its effective medium model become significantly different because the effective medium approach cannot capture the microscopic evanescent and propagating waves in different dielectric layers and tunneling effects [16,17]. These works are focused on the one-dimensional (1D) dielectric multilayers, in which evanescent waves occur only under large incident angles.

In 2D and 3D dielectric composite structures, the scenario is totally different. Inhomogeneity in 2D and 3D models always produces scattering fields, which contain both propagating waves [black arrows in Fig. 1(a)] and evanescent waves [red wavy arrows in Fig. 1(a)] [28]. Such evanescent waves produce rapidly varying evanescent fields nearby the interfaces of dielectric inhomogeneities even in deep-subwavelength structures. Since the evanescent fields occur in a very small area, and simultaneously increase and decrease at different locations, they are averaged out in the traditional Maxwell Garnett EMT.

In this work, we investigate 2D and 3D all-dielectric composite structures at the deep-subwavelength scale. We find that no matter how small the inhomogeneities are (e.g., even at a scale $<\lambda /100$, $\lambda $ is the wavelength in free space), evanescent fields will always occur nearby the interfaces of dielectric inhomogeneities. Because of these inevitable evanescent fields, Maxwell Garnett EMT breaks down when tiny absorptive constituents are placed into the system. We find that the absorption and transmission of the whole structure rely strongly on the positions of the tiny lossy inclusions because they could experience totally different local fields at different positions. Since the traditional Maxwell Garnett EMT averages out the evanescent fields, it fails to predict such position-dependent characteristics. By taking the evanescent waves into consideration, we have developed a correction to the Maxwell Garnett EMT, which can accurately predict the position-dependent transmission and absorption for both 2D and 3D models. Moreover, we predict an interesting phenomenon of “invisible” loss induced by the prohibition of electric fields that appears besides the epsilon-near-zero (ENZ) inclusions. Our work thus reveals a mechanism to control the bulk properties of photonic composite materials without changing the composition. This mechanism is beyond the description of the traditional EMT.

## 2. BREAKDOWN OF MAXWELL GARNETT THEORY AND THE CORRECTION

We first consider the deep-subwavelength 2D model illustrated in Fig. 2(a) under the illumination of transverse-magnetic (TM, out-of-plane magnetic fields) polarized waves. Generally, incident waves will be scattered into various directions by different inclusions. Such scattered waves are usually considered to be propagating waves (black arrows). However, scattered evanescent waves (red wavy arrows) will also emerge if the sizes of the inclusions are comparable to or smaller than the wavelength, as sketched in Fig. 1(a). For visualization, we simulate the wave propagation in a dielectric composite consisting of a dielectric host (relative permittivity ${\epsilon}_{h}=2$, width $w$, height $h=0.8w$) and three different dielectric inclusions by using the finite-element software COMSOL Multiphysics (COMSOL, Inc., Burlington, MA, USA). The three inclusions are of elliptical, circular, and square cross sections, and exhibit relative permittivities of 1, 5, and 3, respectively. Periodic boundary conditions are set on the upper and lower boundaries, and a TM-polarized plane wave with an electric-field amplitude of 1 V/m and a wavelength of $\lambda =125h$ is incident from the free space on the left side. Figures 1(b)–1(d), respectively, present the snapshots of $z$-component ${E}_{z}$, $x$-component ${E}_{x}$, and amplitude $|\mathbf{E}|$ of the electric fields. There are clearly rapidly varying electric fields nearby the inclusions (dielectric interfaces). These evanescent waves are evanescent in the forward and backward directions (i.e., the $z$ direction), but capable of transferring energy flux from high-$\epsilon $ regions to low-$\epsilon $ regions along the perpendicular directions (i.e., the $x$ direction) [28]. They can be understood as a direct consequence of the electromagnetic field boundary conditions. Generally, at this deep-subwavelength scale, the electric fields are enhanced in the low-$\epsilon $ inclusions (the elliptical one), but weakened in the high-$\epsilon $ inclusions (the circular and square ones), as seen in Fig. 1(d).

According to the Maxwell Garnett EMT, we can calculate the effective permittivity ${\epsilon}_{\mathrm{eff}}$ of a $d$-dimensional composite based on [3]

Nevertheless, the Maxwell Garnett EMT will become inaccurate and even break down completely when the composite contains tiny dissipative inclusions, whose sizes are comparable to or smaller than the decay length of the evanescent fields. When such tiny absorptive inclusions are placed close to the large inclusions, they will experience strongly enhanced or weakened fields instead of the averaged fields, which will lead to position-dependent absorption and transmission characteristics. In terms of absorption, the deviation from the Maxwell Garnett EMT could be significantly large, as we will demonstrate below.

For simplicity, we consider an example of a 2D dielectric composite consisting of a dielectric host (relative permittivity ${\epsilon}_{h}$, width $w$, height $h=0.8w$) and an inclusion with circular cross section (relative permittivity ${\epsilon}_{1}$, radius ${r}_{1}=0.15w$), as shown in Fig. 2. A TM-polarized wave with wavelength of $\lambda =125h$ is incident from the free space on the left side. Figures 2(a) and 2(c) present the simulated $|\mathbf{E}|$-distributions for the case with ${\epsilon}_{h}=2$ and ${\epsilon}_{1}=5$, and the case with ${\epsilon}_{h}=5$ and ${\epsilon}_{1}=2$, respectively. In Fig. 2(a), it is seen that the electric fields on the upper/lower side of the inclusion (e.g., position 1) are enhanced, while those on the left/right side are weakened (e.g., positions 2 and 3). In Fig. 2(c), the situation is just the opposite. The existence of such evanescent fields is guaranteed by the electromagnetic field boundary conditions on the interface of the inclusions. Next, we provide an understanding based on the quasi-static model. Since the inclusion is in the deep-subwavelength scale, we assume a background uniform electric field ${\mathbf{E}}_{0}$ (along the $x$ direction) in the host without inclusions, which is induced by the electric field of incidence in simulations. Then, the electric field inside the circular inclusion will be ${\mathbf{E}}_{1}=\frac{2{\epsilon}_{h}}{{\epsilon}_{h}+{\epsilon}_{1}}{\mathbf{E}}_{0}$ [29]. When ${\epsilon}_{1}>{\epsilon}_{h}$ (or ${\epsilon}_{1}<{\epsilon}_{h}$), we have ${\mathbf{E}}_{1}<{\mathbf{E}}_{0}$ (${\mathbf{E}}_{1}>{\mathbf{E}}_{0}$), indicating weakened (or enhanced) electric fields inside the high-$\epsilon $ (or low-$\epsilon $) inclusions, as observed in Figs. 1(d), 2(a), and 2(c). Considering the continuity boundary condition at the host-inclusion interface, we find that the electric fields in the host nearby the interface are

For demonstration, in Fig. 2(a), we add an additional dielectric inclusion with a material loss (relative permittivity ${\epsilon}_{a}=2+i$, radius ${r}_{a}={r}_{1}/6$), as illustrated by dashed circles at positions 1–4. The additional inclusion is much smaller than the original one, so that the original evanescent fields are not largely disturbed and the additional inclusion can experience enhanced or weakened local fields instead of the averaged fields. The absorptance by the composite is plotted in Fig. 2(b) as a function of the working wavelength when the additional inclusion successively moves from position 1 to position 4. We can see that the wave absorption for the additional inclusion at position 1 is much larger than that at other positions, because the electric field induced by evanescent scattering waves is largest at position 1. This clearly shows the dependence of wave absorption on the positions of the tiny inclusion. However, the traditional EMT [i.e., Eq. (1)] ignores the evanescent fields and thus gives the same absorptance [black solid lines in Fig. 2(b)], which almost coincides with the result for the case of position 4 where evanescent waves almost disappear. The traditional EMT cannot capture the details of evanescent fields at the deep-subwavelength scale and thus fails in correctly describing such a position-dependent absorption.

Interestingly, by taking the evanescent fields into consideration, we can modify the formula of the traditional EMT, to give correct description of the position-dependent absorption and transmission characteristics. The effective permittivity of the $d$-dimensional composite containing $M$ large inclusions and ${M}^{\prime}$ tiny inclusions can be calculated based on the corrected formula as

From the field-distribution in Fig. 2(a), we find that the correction factor $\beta $ at positions 1–4 is around 1.386, 0.7158, 0.7165, and 1.028, respectively. The $\beta $ at positions 2 and 3 is nearly the same, because the composite lies in the deep subwavelength scale. If the working wavelength tends to be infinitely long, we will get exactly the same $\beta $ under the electrostatic limit [see Eq. (2)]. We also see that the $\beta $ at position 4 is near unity, indicating that the evanescent fields are mostly localized in very limited areas. Since Eq. (3) becomes Eq. (1) under $\beta =1$, the absorption of the case at position 4 is close to the absorption predicted by the traditional EMT, as seen in Fig. 2(b). According to the corrected EMT [Eq. (3)], we plot the absorptance for the cases with the tiny inclusion at positions 1–3, as shown by the dashed lines in Fig. 2(b). We see that the results coincide with the simulation results quite well, demonstrating the validity of the proposed correction in the EMT.

We also reanalyze the composite with ${\epsilon}_{h}=5$ and ${\epsilon}_{1}=2$ in Fig. 2(c) by adding an additional dielectric inclusion with a material loss (${\epsilon}_{a}=5+i$, ${r}_{a}={r}_{1}/6$). The correction factor $\beta $ at positions 1–4 is found to be around 0.6486, 1.254, 1.257, and 0.9730, respectively. The absorptance of the composite as a function of the working wavelength is calculated based on the simulations of the actual composite, traditional EMT [Eq. (1)], and corrected EMT [Eq. (3)], as shown by the dots, solid lines, and dashed lines in Fig. 2(d), respectively. We see that the traditional EMT fails to correctly describe the wave absorption when the additional tiny inclusion is close to the original large inclusion. Interestingly, the proposed corrected EMT can give an accurate description. We note that the deviation for the cases at positions 2 and 3 is mainly caused by the inhomogeneity of fields where the additional inclusion lies. If we further reduce the size of the additional inclusion and increase the wavelength, the deviation will become smaller.

The proposed correction of the EMT can be applied to dielectric composites containing arbitrarily shaped inclusions. As an example, we consider a dielectric composite consisting of a dielectric host (${\epsilon}_{h}=2$, width $w$, height $h=0.8w$) embedded with two cloud-like inclusions (${\epsilon}_{1}=5$, ${\epsilon}_{2}=1$) and three additional tiny inclusions (having the same relative permittivity ${\epsilon}_{a}=2+0.5i$, radius ${r}_{a}=0.015w$) placed at positions 1–3, as illustrated in Fig. 3(a). Figure 3(b) presents the simulated $|\mathbf{E}|$-distribution in the absence of the additional inclusions under the illumination of a TM-polarized wave with $\lambda =125h$, showing rapidly varying fields nearby the two cloud-like inclusions. Based on this field-distribution, the correction factors at positions 1–3 are evaluated as 1.246, 1.308, and 0.697, respectively. The absorptance of this composite is calculated as the function of wavelength based on the simulations of the actual composite, traditional EMT [Eq. (1)], and corrected EMT [Eq. (3)], as shown by the dots, solid lines, and dashed lines in Fig. 3(c), respectively. The results clearly demonstrate the breakdown of the traditional EMT, and the validity of the corrected EMT in such a complex composite. In addition, in Fig. 3(d), we compare the transmittance through the effective media based on the traditional EMT (black lines) and the corrected EMT (red lines) by assuming $N$ layers of unit cells along the propagation direction. The results clearly show that the significant difference in absorption could lead to a large deviation in the transmission for large samples. Here, we note that because the thickness of the composite sample ($\sim 100\lambda $) is much larger than the working wavelength, the interference of multiple reflections between the interfaces is omitted in the calculation.

## 3. INVISIBLE LOSS INDUCED BY EVANESCENT FIELDS

The position-dependent characteristic can lead to an interesting phenomenon (i.e., the disappearance of absorption in dielectric composites with absorptive constituents). This is impossible from the viewpoint of traditional EMT, because all inclusions contribute to the effective parameters. Interestingly, we find that when tiny absorptive inclusions are placed where $\beta \to 0$, the waves cannot “see” them. As a result, the effective permittivity of the whole composite will be totally independent of these tiny absorptive inclusions [see Eq. (3)], thus leading to the disappearance of total absorption. Equation (2) indicates that the area with $\beta \to 0$ can actually be obtained if the composite contains ENZ inclusions, which is also denoted as “side scattering shadows” [30].

As an example, we consider a 2D composite consisting of a dielectric host (${\epsilon}_{h}=2$, width $w$, height $h=0.8w$) embedded with a rectangular ENZ inclusion (${\epsilon}_{1}=0.001$, width $0.7w$, height $0.1w$) and two slabs of lossy dielectric (${\epsilon}_{a}=2+i$, width $0.35w$, height $0.01w$) coated on the upper and lower surfaces of the ENZ inclusion, as sketched in Fig. 4(a). Figure 4(b) displays the simulated $|\mathbf{E}|$-distribution in the absence of the two lossy inclusions when a TM-polarized wave with $\lambda =125h$ is normally incident from the free space on the left side. Clearly, due to the matching of the displacement fields inside and outside the ENZ medium, the electric fields become extremely small on the upper and lower sides of the ENZ inclusion. From the average field in the area where the lossy inclusions lie, we evaluate the correction factor $\beta $ as 0.151, which is significantly smaller than unity. Figure 4(c) presents the absorptance of this composite as a function of the working wavelength. As expected, the wave absorption in the actual composite (dots) is negligibly small, which is well predicted by the corrected EMT (dashed lines). The absorption predicted by the traditional EMT (solid lines) shows an obvious deviation from the actual absorption. Moreover, we plot the transmittance through the effective media based on the traditional EMT (black lines) and the corrected EMT (red lines) with $N$ layers of unit cells in Fig. 4(d). We can see that the transmission in both cases oscillates over $N$ as the result of Fabry–Perot resonances. The transmission predicted by the traditional EMT (black lines) decreases quickly when $N$ increases. At the same time, the transmission oscillation predicted by the corrected EMT remains almost unchanged (red lines), indicating the disappearance of loss.

## 4. 3D MODELS

Besides 2D models, the breakdown of traditional Maxwell Garnett EMT also applies to 3D models. Figure 5(a) shows a practical 3D model consisting of a sphere of silicon (Si, ${\epsilon}_{1}=12$, radius 15 nm) and eight spheres of indium tin oxide (ITO, radius 3 nm) in a host of silica (${\mathrm{SiO}}_{2}$, ${\epsilon}_{h}=2.1$, side length 50 nm). The eight ITO spheres are evenly distributed along the equator line of the Si sphere. The relative permittivity of the ITO sphere is $\epsilon (\omega )={\epsilon}_{\infty}-\frac{{\omega}_{p}^{2}}{\omega (\omega +i\mathrm{\Gamma})}$ [32], where ${\epsilon}_{\infty}=3\text{.}94$, ${\omega}_{p}=\sqrt{5.97}\times {10}^{15}\text{\hspace{0.17em}}\mathrm{Hz}$, $\mathrm{\Gamma}=1.88\times {10}^{14}\text{\hspace{0.17em}}\mathrm{Hz}$, and $\omega $ is the angular frequency. A plane wave with electric fields polarized along the $x$ direction is normally incident from the free space on the left side. Figure 5(b) presents the simulated $|\mathbf{E}|$-distribution in the absence of the ITO spheres when the working wavelength is 1400 nm, showing enhanced fields nearby the poles and weakened fields nearby the equator line of the Si sphere. Similar to Eq. (2) for the 2D model, we can also evaluate the electric fields in this 3D model under the electrostatic limit as [29]

Since ${\epsilon}_{1}>{\epsilon}_{h}$ in the model in Fig. 5(a), we have ${\mathbf{E}}_{h}^{\text{equator}}<{\mathbf{E}}_{0}$ and ${\mathbf{E}}_{h}^{\mathrm{pole}}>{\mathbf{E}}_{0}$, as observed in Fig. 5(b). This indicates the rapid variation of electric fields nearby the Si sphere (i.e., the existence of evanescent fields). Likewise, if additional tiny inclusions experiencing the evanescent local fields exit, the traditional Maxwell Garnett EMT will fail to describe such location-dependent situations.

For verification, we have calculated the absorptance as a function of the working wavelength, as shown in Fig. 5(c). The dots and solid lines denote the results obtained through the simulation of the actual composite and the traditional EMT [Eq. (1)], respectively. Clearly, the effective medium model overestimates the absorption because the fields experienced by the ITO spheres are weakened. Interestingly, the proposed correction of EMT in Eq. (3) can accurately describe the situation for such 3D models. Based on the field distribution, we find out the correction factor is 0.66. Then, we calculate the absorptance via the corrected EMT according to Eq. (3), and plot the results as dashed lines in Fig. 5(c), showing a good match with the actual composite simulation results. Moreover, we plot the transmittance on a log scale [i.e., $\mathrm{log}(T)$] through the effective media based on the traditional EMT (black lines) and the corrected EMT (red lines) with $N$ layers of unit cells at the wavelength of 1400 nm, as displayed in Fig. 5(d). The interference of multiple reflections is omitted because the total thickness ($\sim 300\text{\hspace{0.17em}}\mathrm{\mu m}$) is much larger than the working wavelength. It is seen that the transmittance predicted by the traditional EMT decreases much faster than that predicted by the corrected EMT.

## 5. DISCUSSION AND CONCLUSION

It is noteworthy that the physical mechanism of the breakdown of EMT in 1D dielectric multilayers [16–27] and 2D/3D dielectric composite structures studied here is fundamentally different. In 1D dielectric multilayers, the EMT breaks down close to the total internal reflection angle originating from tunneling effects of evanescent waves. Under the critical incident angle, the waves become evanescent in low-$\epsilon $ layers, but remain propagating in high-$\epsilon $ layers. Since the layers are deep-subwavelength, the incident waves may still propagate through the multilayer via tunneling, whereas the EMT does not capture this physics, thus leading to the failure of the EMT [16,17]. Nevertheless, in our proposed 2D/3D structures, the fundamental origin for the breakdown of the EMT is the dramatically varying evanescent fields induced by the field-matching condition on the surfaces of the inhomogeneities. When there are tiny absorptive inclusions experiencing such varying local fields, the macroscopic properties (e.g., reflection, transmission, and absorption) of the composite structures become sensitive to the positions of those tiny inclusions. Such varying local fields can be comprehended as the excitation of high-order modes [33] instead of a dipole–dipole interaction [34]. However, in the traditional EMT description, these details are averaged out, thus leading to the breakdown of the EMT. We note that the breakdown in the 2D/3D structures does not rely on the angle of incidence, which can be observed even under normal incidence as shown above.

The continuity of the electric displacement at the inclusion–host interfaces plays an important role in generating the dramatically varying evanescent fields at this deep-subwavelength scale. For instance, in 3D models, the electric fields inside inclusions are nearly uniform. Since the inclusion–host interface is spherical, it is parallel to the fields at some places, but perpendicular to them at other places. This leads to rapidly varying fields in the host nearby the inclusion–host interface. This scenario can also be seen in 2D models for the TM polarization, as we have demonstrated above. However, we also note that for transverse-electric (TE, out-of-plane electric fields) polarization in 2D models, the electric fields nearby the inclusion–host interface are always parallel to the interface. In this case, the traditional Maxwell Garnett EMT is still valid. Additional details can be found in the supplementary materials of Ref. [31].

Alternatively, the breakdown of the EMT can also be understood from the mode interactions. In the deep-subwavelength scale, the dipole mode of particles dominates, while high-order modes are generally negligibly small. For subwavelength objects far from each other, the dipole approximation is valid, which is the basis of the traditional EMT [2]. For objects close to each other, however, the contribution from high-order modes would be dramatic [33], thus leading to the failure of the traditional EMT.

We note that the discovered position-dependent transmission/absorption characteristics at the deep-subwavelength scale are beyond the extended EMT [35–37]. The extended EMT is usually used to deal with the composite with components rather than that subwavelength, where high-order terms must be considered. However, in our deep-subwavelength structures, the high-order terms would be negligibly small. In such deep-subwavelength models, the extended EMT will be consistent with the traditional EMT.

Finally, it is also worth noting that for simple structures like cylinders and spheres, the correction factor $\beta $ can be analytically evaluated based on Eqs. (2) and (4) without using numerical simulations, which leads to $\beta $ in the range of 0–2 (or 0–1.5 or 3) for the 2D cylindrical (or 3D spherical) models. Based on this range of values, we can immediately obtain the range of potential deviation for the EMT estimation. This information is valuable in many situations. For example, in complicated structures (e.g., Figs. 3 and 4), although numerical analysis is needed to precisely compute $\beta $ because there are no longer simple analytical solutions, our theory can still serve as a guide to design and control absorption without changing the composition of the material. This is unimaginable from the viewpoint of traditional EMT, which is deeply believed in optics. We believe, however, that the numerical simulation will not ruin the value and universality of our proposed corrected EMT.

In summary, we have considered the model of dielectric composites with absorptive constituents. Such a description is generally valid in many circumstances where the absorption is mainly induced by some tiny particles or molecules in the system. Because the region of absorption is determined by the positions of the absorptive constituents, such a case can maximize the difference induced by the evanescent fields at the deep-subwavelength scale. The breakdown of the traditional EMT is inevitable because it simply averages out the evanescent fields and ignores their feature. A correction by taking the distribution of evanescent fields into consideration can significantly increase the accuracy of the EMT prediction. Although traditional wisdom tells us that dielectric structures at the deep-subwavelength scale can be well predicted by the EMT based on homogenized fields, our finding reveals an intriguing exceptional case where the traditional EMT fundamentally breaks down. We have demonstrated that a microscopic variation of dielectric structure at the deep-subwavelength scale can also lead to a dramatic difference in bulk behaviors, even when the composition of the composite is fixed. Actually, the proposed configuration is common in biomedicine and molecular biology. Nanoparticles are usually used to defect biomolecules and diagnostic assay, forming a configuration consisting of large nanoparticles and nearby absorptive biomolecules [38,39]. We therefore believe that in addition to the importance of understanding the EMT, our work will also be very interesting and will have an effect on advanced photonics, bioscience, and practical applications.

## Funding

National Key Research and Development Program of China (2017YFA0303702); National Natural Science Foundation of China (11704271, 11974176, 61671314); Natural Science Foundation of Jiangsu Province (BK20170326); Priority Academic Program Development of Jiangsu Higher Education Institutions.

## Disclosures

The authors declare no conflicts of interest.

## REFERENCES

**1. **J. C. M. Garnett, “Colours in metal glasses and in metallic films,” Philos. Trans. R. Soc. A **203**, 385–420 (1904). [CrossRef]

**2. **V. A. Markel, “Introduction to the Maxwell Garnett approximation: tutorial,” J. Opt. Soc. Am. A **33**, 1244–1256 (2016). [CrossRef]

**3. **T. C. Choy, *Effective Medium Theory* (Oxford University, 1999).

**4. **K. Dolgaleva and R. W. Boyd, “Local-field effects in nanostructured photonic materials,” Adv. Opt. Photon. **4**, 1–77 (2012). [CrossRef]

**5. **W. Cai and V. Shalaev, *Optical Metamaterials: Fundamentals and Applications* (Springer, 2009).

**6. **P. A. Belov, R. Marqués, S. I. Maslovski, I. S. Nefedov, M. Silveirinha, C. R. Simovski, and S. A. Tretyakov, “Strong spatial dispersion in wire media in the very large wavelength limit,” Phys. Rev. B **67**, 113103 (2003). [CrossRef]

**7. **L. Shen, T. Yang, and Y. Chau, “50/50 beam splitter using a one-dimensional metal photonic crystal with parabolalike dispersion,” Appl. Phys. Lett. **90**, 251909 (2007). [CrossRef]

**8. **R. Pollard, A. Murphy, W. Hendren, P. Evans, R. Atkinson, G. Wurtz, A. Zayats, and V. Podolskiy, “Optical nonlocalities and additional waves in epsilon-near-zero metamaterials,” Phys. Rev. Lett. **102**, 127405 (2009). [CrossRef]

**9. **C. R. Simovski, P. A. Belov, A. V. Atrashchenko, and Y. S. Kivshar, “Wire metamaterials: Physics and applications,” Adv. Mater. **24**, 4229–4248 (2012). [CrossRef]

**10. **J. Luo, H. Chen, B. Hou, P. Xu, and Y. Lai, “Nonlocality-induced negative refraction and subwavelength imaging by parabolic dispersions in metal-dielectric multilayered structures with effective zero permittivity,” Plasmonics **8**, 1095–1099 (2013). [CrossRef]

**11. **R. L. Chern, “Spatial dispersion and nonlocal effective permittivity for periodic layered metamaterials,” Opt. Express **21**, 16514–16527 (2013). [CrossRef]

**12. **M. G. Silveirinha, “Nonlocal homogenization model for a periodic array of ϵ-negative rods,” Phys. Rev. E **73**, 046612 (2006). [CrossRef]

**13. **A. Alù, “First-principles homogenization theory for periodic metamaterials,” Phys. Rev. B **84**, 075153 (2011). [CrossRef]

**14. **J. Luo, Y. Yang, Z. Yao, W. Lu, B. Hou, Z. H. Hang, C. T. Chan, and Y. Lai, “Ultratransparent media and transformation optics with shifted spatial dispersions,” Phys. Rev. Lett. **117**, 223901 (2016). [CrossRef]

**15. **S. Li, Y. Wang, W. Zhang, W. Lu, B. Hou, J. Luo, and Y. Lai, “Observation of wide-angle impedance matching in terahertz photonic crystals,” New J. Phys. **22**, 023033 (2020). [CrossRef]

**16. **H. H. Sheinfux, I. Kaminer, Y. Plotnik, G. Bartal, and M. Segev, “Subwavelength multilayer dielectrics: ultrasensitive transmission and breakdown of effective-medium theory,” Phys. Rev. Lett. **113**, 243901 (2014). [CrossRef]

**17. **S. V. Zhukovsky, A. Andryieuski, O. Takayama, E. Shkondin, R. Malureanu, F. Jensen, and A. V. Lavrinenko, “Experimental demonstration of effective medium approximation breakdown in deeply subwavelength all-dielectric multilayers,” Phys. Rev. Lett. **115**, 177402 (2015). [CrossRef]

**18. **A. Andryieuski, A. V. Lavrinenko, and S. V. Zhukovsky, “Anomalous effective medium approximation breakdown in deeply subwavelength all-dielectric photonic multilayers,” Nanotechnology **26**, 184001 (2015). [CrossRef]

**19. **H. H. Sheinfux, I. Kaminer, A. Z. Genack, and M. Segev, “Interplay between evanescence and disorder in deep subwavelength photonic structures,” Nat. Commun. **7**, 12927 (2016). [CrossRef]

**20. **M. Coppolaro, G. Castaldi, and V. Galdi, “Effects of deterministic disorder at deeply subwavelength scales in multilayered dielectric metamaterials,” Opt. Express **28**, 10199–10209 (2020). [CrossRef]

**21. **D. V. Novitsky, A. S. Shalin, and A. Novitsky, “Nonlocal homogenization of PT-symmetric multilayered structures,” Phys. Rev. A **99**, 043812 (2019). [CrossRef]

**22. **M. A. Gorlach and M. Lapine, “Boundary conditions for the effective-medium description of subwavelength multilayered structures,” Phys. Rev. B **101**, 075127 (2020). [CrossRef]

**23. **X. Lei, L. Mao, Y. Lu, and P. Wang, “Revisiting the effective medium approximation in all-dielectric subwavelength multilayers: breakdown and rebuilding,” Phys. Rev. B **96**, 035439 (2017). [CrossRef]

**24. **V. Popov, A. V. Lavrinenko, and A. Novitsky, “Operator approach to effective medium theory to overcome a breakdown of Maxwell Garnett approximation,” Phys. Rev. B **94**, 085428 (2016). [CrossRef]

**25. **A. Maurel and J. Marigo, “Sensitivity of a dielectric layered structure on a scale below the periodicity: a local homogenized model,” Phys. Rev. B **98**, 024306 (2018). [CrossRef]

**26. **M. Coppolaro, G. Castaldi, and V. Galdi, “Anomalous light transport induced by deeply subwavelength quasiperiodicity in multilayered dielectric metamaterials,” Phys. Rev. B **102**, 075107 (2020). [CrossRef]

**27. **H. H. Sheinfux, Y. Lumer, G. Ankonina, A. Z. Genack, G. Bartal, and M. Segev, “Observation of Anderson localization in disordered nanophotonic structures,” Science **356**, 953–956 (2017). [CrossRef]

**28. **J. Luo, W. Lu, Z. Hang, H. Chen, B. Hou, Y. Lai, and C. T. Chan, “Arbitrary control of electromagnetic flux in inhomogeneous anisotropic media with near-zero index,” Phys. Rev. Lett. **112**, 073903 (2014). [CrossRef]

**29. **D. J. Griffiths, *Introduction to Electrodynamics*, 3rd ed. (Prentice Hall, 1999).

**30. **J. Song, J. Luo, and Y. Lai, “Side scattering shadow and energy concentration effects of epsilon-near-zero media,” Opt. Lett. **43**, 1738–1741 (2018). [CrossRef]

**31. **T. Dong, J. Luo, H. Chu, X. Xiong, and Y. Lai, “Breakdown of Maxwell Garnett theory due to evanescent fields at deep-subwavelength scale,” arXiv:2012.10088 (2020).

**32. **A. Capretti, Y. Wang, N. Engheta, and L. D. Negro, “Enhanced third-harmonic generation in Si-compatible epsilon-near-zero indium tin oxide nanolayers,” Opt. Lett. **40**, 1500–1503 (2015). [CrossRef]

**33. **J. B. Pendry, A. I. Fernández-Domínguez, Y. Luo, and R. Zhao, “Capturing photons with transformation optics,” Nat. Phys. **9**, 518–522 (2013). [CrossRef]

**34. **B. X. Wang and C. Y. Zhao, “Near-resonant light transmission in two-dimensional dense cold atomic media with short-range positional correlations,” J. Opt. Soc. Am. B **37**, 1757–1768 (2020). [CrossRef]

**35. **R. Ruppin, “Evaluation of extended Maxwell-Garnett theories,” Opt. Commun. **182**, 273–279 (2000). [CrossRef]

**36. **S. Chui and L. Hu, “Theoretical investigation on the possibility of preparing left-handed materials in metallic magnetic granular composites,” Phys. Rev. B **65**, 144407 (2002). [CrossRef]

**37. **Y. Wu, J. Li, Z. Q. Zhang, and C. T. Chan, “Effective medium theory for magnetodielectric composites: beyond the long-wavelength limit,” Phys. Rev. B **74**, 085111 (2006). [CrossRef]

**38. **N. C. Tansil and Z. Gao, “Nanoparticles in biomolecular detection,” Nano Today **1**, 28–37 (2006). [CrossRef]

**39. **A. van Reenen, A. M. de Jong, J. M. den Toonder, and M. W. Prins, “Integrated lab-on-chip biosensing systems based on magnetic particle actuation–a comprehensive review,” Lab Chip **14**, 1966–1986 (2014). [CrossRef]