## Abstract

Metal nanorod arrays exhibit hyperbolic dispersion and optical nonlocality under certain conditions. Therefore, their optical behaviors can hardly be expressed by incident-angle-independent effective permittivity. Here we extract effective permittivity of silver nanorod arrays with diameters of 4 nm, 12 nm, and 20 nm by polarized transmission method in the visible range. The incident angles are chosen from 20° to 60° to study the influence of optical nonlocality on permittivity. We demonstrate how the diameter of the nanorods can control the effective permittivity beyond the effective medium theory. The results suggest that the effective permittivity gradually loses its accuracy as the diameter increases due to the optical nonlocality. Our experiment verifies that ultrathin nanorod arrays can resist the fluctuations caused by changes in incident angle. We also extract *k*-dependent effective permittivity of nanorods with larger diameters.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Hyperbolic metamaterials (HMMs) are anisotropic materials that one of the in-plane and out-of-plane permittivity components is positive, and the other is negative [1]. They reveal plenty of unique optical effects for their extreme anisotropy and epsilon-near-zero (ENZ), such as enhanced nonlinearity [2–4], perfect absorption [5], engineered thermal radiation [6], and magneto-optical effects [7]. HMMs are usually realized by metal nanorod arrays or alternating layers of dielectrics and metals [1]. Generally, their effective permittivity is quantitatively obtained via the Maxwell-Garnet effective medium theory (EMT) in recent reports [8–11]. This approximation has been proved to be valid in many cases [12,13]. However, local EMT neglects the optical nonlocality. Therefore, the obtained effective permittivity has obvious errors in predicting some special phenomena of HMMs, such as anomalous transmittance near the critical angle of total internal reflection [14,15], angle-dependent photonic gaps in layered HMMs, and surface wave beyond the EMT [16,17].

Numerous studies have been carried out around the measurement of the effective permittivity of HMMs in recent years. The measurement methods include reflectivity method, S parameter method, and ellipsometry method [18–21]. However, the accuracy of these methods is often limited since it is difficult to model the permittivity of HMMs. T. Tumkur et al. modeled the permittivity of alternating layers of dielectric and metal by hybrid Gaussian-Lorentz oscillators, then extracted it from ellipsometry parameters. This method obtained both in-plane and out-of-plane permittivity. However, the results are noisy and do not satisfy the Kramers-Kronig relationship [19]. O. Yermakov et al. modeled the in-plane conductivity of elliptical gold nanodisk arrays by a dipole resonance model, then extracted it from S parameter measured in the far-field. This method obtained accurate in-plane conductivity of type-I HMM (Re(*ɛ*_{||}) > 0 and Re(*ɛ*_{⊥}) < 0). The out-of-plane conductivity was not achieved by this method [21].

Recently, Cheng Zhang et al. modeled the out-of-plane permittivity of alternating layers of Ag and Ta_{2}O_{5} by B-spline free parameters and modeled the in-plane permittivity by oscillators. Then they extracted the permittivity from ellipsometry parameters measured through total internal reflection at 60°. The method was proved to be robust [22]. James Dilts et al. modeled the out-of-plane permittivity of three kinds of type-II HMMs (Re(*ɛ*_{||}) < 0 and Re(*ɛ*_{⊥}) > 0) by B-spline free parameters and calculated the in-plane permittivity by the EMT. They then obtained the out-of-plane permittivity from ellipsometry parameters measured at 65°, 70°, and 75°. This method was effective in reducing the mean square error [23]. However, these methods were performed over a small range of incident angles, so the applicability at other incident angles may be limited due to the optical nonlocality. Changes in the incident angle may lead to a set of different effective permittivity. As a result, the error of the effective permittivity extracted from a single incident angle is obviously lower than that obtained from multiple incident angles [22]. Recent reports have also demonstrated that the optical behavior of metamaterials cannot be explained by incident-angle-independent effective permittivity, even in simple structures [24]. On account of the existence of optical nonlocality, the permittivity of HMMs is related to the incident angle [25–28]. Optical nonlocality plays an essential role in describing the optical behaviors of HMMs, especially near the ENZ position for the enhancement of local electric fields. These optical characteristics can hardly be expressed by local approximations [29–33], such as the beam splitting near the Dirac point of alternating layers of dielectric and metal [34]. More in-depth research on optical nonlocality can help achieving flexible manipulation of optical spectra and electric fields [35,36].

Metal nanorod arrays are a kind of type-I HMMs. The special electric field distribution in nanorod arrays results in strong optical nonlocality, which makes it difficult to characterize the material at multiple incident angles [37,38]. However, metal nanorod arrays have potential applications in super-resolution and imaging [39]. Accurately predicting the transmission spectrum is beneficial. To obtain accurate optical permittivity of nanorod arrays, and find the influence of nonlocality on the nanorods of different sizes, in this work, we extract the effective permittivity of three nanorod arrays. Their diameters are 4 nm, 12 nm, and 20 nm, respectively. They possess the same volume fraction of silver. We model the in-plane and out-of-plane permittivity of metal nanorod arrays by B-spline free parameter. Then we obtain the complex permittivity from the transmission spectra for s and p incident light measured from 20° to 60°. Compared with [22] and [23], which used the same method to extract the effective permittivity of HMMs, we extend this method to type-I HMMs for the first time. And we perform this method over a wider range of incident angles to enhance the influence of optical nonlocality. We also use this method to extract *k*-dependent permittivity, which is more reasonable for expressing the optical characteristics of HMMs. We find the accuracy of the *k*-independent effective permittivity is significantly improved in the nanorod arrays of small diameter. Our experiment confirms the applicability of the effective permittivity to the ultrathin nanorod array. For the nanorod arrays of large diameters, we extract *k*-dependent effective permittivity to express their transmittance.

## 2. Representation of silver nanorod arrays

An ultrathin silver nanorod array, *d* = 4 nm in diameter, *S* = 9 nm in space, and *L* = 120 nm in length, surrounded by aluminum oxide is chosen to be analyzed in this work. The schematic diagram of the sample is illustrated in Fig. 1(a). Two other samples with *d* = 12 nm, *S* = 27 nm and *d* = 20 nm, *S* = 45 nm are also analyzed for comparison. The volume fraction of silver in the three samples is *ρ *= π*d*^{2}/(4*S*^{2}) = 15.51%. The transmission spectrum of the metal nanorod arrays is dominated by two distinct absorption peaks for p incident light. And only the absorption peak at shorter wavelengths appears for s polarization incident light [40,41]. Transmission spectra of the three samples are calculated by finite element method (FEM) and presented in Fig. 1(b). The two absorption peaks are located at around 400 nm and 600 nm respectively in the transmission spectrum.

To quantitatively characterize the optical behavior of metal nanorod arrays, the nanostructures are usually treated as homogeneous, uniaxial films. The optical responses of the film can be described by a set of orthogonal complex permittivity *ɛ*_{⊥} and *ɛ*_{||}, where *ɛ*_{⊥} represents the permittivity tensor perpendicular to the surface of the sample, and *ɛ*_{||} represents the permittivity tensor parallel to the surface of the sample. The permittivity tensor is expressed as follows:

*ɛ*

_{⊥}and

*ɛ*

_{||}are usually calculated by the EMT method. Two kinds of EMTs for metal nanorod arrays are considered in this study. A commonly used local EMT is shown in Eqs. (2 and 3) [38,42]. No nonlocal element is included in this method.

*ρ*is the volume fraction of silver,

*ɛ*

_{Ag}and

*ɛ*

_{Al2O3}are the complex permittivity of silver and Al

_{2}O

_{3}, respectively. In our work,

*ɛ*

_{Ag}is derived from [43] and

*ɛ*

_{Al2O3}from [44].

*λ*is the wavelength of incident light. We refer to this method as size-independent EMT for brevity.

Another kind of EMT related to the specific size of the material is also calculated for comparison. We refer to it as the size-dependent EMT. Partly similar to the nonlocal models, the size-dependent EMT considers the interaction between nanoparticles roughly. However, the influence of the incident angle is neglected. The size-dependent EMT can be shown as follows [45]

*σ*(= || or ⊥) denotes the direction.

*α*is the polarizability of the nanorod arrays. γ

_{σ}*is the tensor indicating the interaction between particles. γ*

_{σ}*varies with the structure. Here γ*

_{σ}_{||}=0.67 and γ

_{⊥}varies with the diameter.

*L*

_{⊥}and

*L*

_{||}are depolarizations related to the shape.

*e*=(1-

*d*

^{2}/

*L*

^{2}) is the particle aspect ratio.

*ɛ*_{||} and *ɛ*_{⊥} calculated by the size-independent EMT are illustrated in Figs. 1(c), (d). The samples exhibit type-I hyperbolic dispersion from 608 nm to the near-infrared region. The three samples have the same volume fraction. Therefore, *ɛ*_{||} and *ɛ*_{⊥} of the three samples are the same when predicted by size-independent EMT. Even so, as shown in Fig. 1(b), transmittance spectra for p polarization incident light of the three samples exhibit considerable difference. One source of the phenomenon is that the size-independent EMT ignores the optical nonlocality of nanorod arrays. The effective permittivity predicted by the size-dependent EMT is also shown in Figs. 1(c) and (d) to show the effect of interaction between the nanorods. The results indicate that the diameter of the nanorods can show an influence on the effective permittivity of the sample. Owing to the effect of interaction between nanorods, the ENZ position shifts toward short wavelengths. The imaginary part of *ɛ*_{⊥} is almost unchanged at 300 nm to 350 nm, and increases with the diameter obviously at wavelengths longer than 350 nm.

Optical nonlocality exists in several periodic metal nanostructures. Their optical parameters show angle sensitivity [24,28,46]. Many efforts have been made to calculate the incident-angle-dependent effective permittivity. For example, the permittivity of a material with optical nonlocality was extended as a function of the wavevector in [46] and [47] as follows:

Some works tried to develop the local EMT for *k*-dependent nonlocal versions [48,49]. For plasmonic nanorod arrays, the nonlocal effective permittivity function is shown as [50].

*q*=

*CρS*4π, and

^{2}/*C*is the Euler constant.

Optical nonlocality is usually related to the special electric field distribution in samples. Previous reports have proved that optical nonlocality in nanorods mainly originates from interference between the electromagnetic waves supported by nanorods [38]. Ordered metal nanorod arrays support an additional electromagnetic wave near the ENZ position. Interference occurs between the main wave and the additional wave, leading to a change in the absorption characteristics and a relatively weak modulation of the effective permittivity. The absorption characteristics of metal nanostructures are closely related to electromagnetic field enhancement. The field distributions of the three kinds of nanorod arrays are illustrated in Fig. 2. Obviously, the electric field distribution changes drastically with the increase in the incident angle. This leads to changes in the absorption characteristics of the samples. Besides that, in nanorods of different sizes, the change in the electric field distribution with the incident angle are also different. When the incident angle increases from 20° to 60°, the electric field becomes more concentrated in nanorods of *d* = 4 nm. However, as shown in Fig. 2(c), another electric field concentration area gradually appears as the incident angle increases to 60° in nanorods of *d* = 20 nm. The special electric field distribution demonstrates that the optical behavior of the nanorod arrays may have angle sensitivity, and the influence of the incident angle on the absorption of nanorod arrays of different sizes is different. Therefore, wavevector-dependent effective optical parameters may be more reasonable when characterizing the optical behavior of nanorod arrays.

## 3. Extraction and discussion effective permittivity

Our purpose is to obtain a set of *ɛ*_{||} and *ɛ*_{⊥}, which can express the transmission spectra for s and p polarization incident light over a wide range of incident angles. Here we regard the silver nanorod arrays as homogeneous, uniaxial films. The thickness of the film is the same as the length of the nanorod arrays. During the calculation, in the first step, we extract *ɛ*_{||} from transmittance spectra for s polarization incident light. In order to verify the applicability of the effective permittivity over a wide range of incident angles, the incident angles are chosen as 20°, 25°, 30°, and 60°. We yield the imaginary part of *ɛ*_{||} and *ɛ*_{⊥} by B-spline free parameter. Their real parts are then solved by the Kramers-Kronig relationship [51]. The initial guesses of the imaginary parts of *ɛ*_{||} and *ɛ*_{⊥} derive from the local size-independent EMT (Eq. (2 and 3) for convenience. Then the transmittance spectra of the homogenized film are calculated by FEM. The mean squared deviation of transmittance (MSDT) is used to compare the difference between calculated data and homogenized models. The MSDT is defined in Eq. (5).

*N*is the number of points used in the calculation.

*T*

_{calculation}and

*T*

_{homogenize}are transmittance spectra obtained from the calculation and homogenized model, respectively.

*λ*is the wavelength of the incident light.

*θ*is the incident angle. Here the incident angles we chose are 20°, 25° 30°, and 60°.

_{i}Regression analysis (Levenberg–Marquardt algorithm) is then used to find the complex permittivity that satisfies the transmittance spectra for s polarization incident light [52]. The Levenberg–Marquardt algorithm is widely used in solving nonlinear least-squares problems. This method has a high calculation speed. Moreover, it can avoid the interference of local minima, which is effective in reducing MSDT. The imaginary part and the real part of *ɛ*_{||} are modulated to minimize the MSDT. Then *ɛ*_{⊥} is solved based on the calculated *ɛ*_{||} to satisfy the transmittance spectra for p polarization incident light. The calculation result is a set of *ɛ*(*λ*). The permittivity is a function of the wavelength.

The best-matched transmittance spectra are illustrated in Fig. 3. The calculation range is selected as 300 nm to 750 nm to cover the two absorption peaks completely. Our purpose here is to describe the optical properties of the material in the strong absorption band accurately, so other wavelengths with weak absorption are not considered. To verify the applicable conditions of the effective permittivity, the MSDT at every incident angle and every polarization state is calculated separately. The values are shown in Table 1. The MSDTs calculated by the size-independent EMT and the size-dependent EMT are shown in Table 2 for comparison. The effective permittivity is significantly more accurate in predicting the transmittance spectra. We set the criteria for judging that the permittivity is reasonable as the MSDT ≤ 2% at all incident angles. It can be observed that the model is in good agreement with the transmittance spectra for s polarization incident light. MSDTs for the three nanorods of different sizes are standards-compliant at all incident angles. However, obvious errors occur when the effective permittivity is used to predict the transmittance spectra for p polarization incident light, especially at wavelengths near the ENZ position. That is, the effective permittivity cannot accurately represent the absorption peaks near 600 nm at all incident angles at the same time. As proved in previous work, the optical nonlocality of this position is the strongest in the entire visible-light band [38]. Therefore, a wavevector-independent effective permittivity cannot express the optical behavior of the material in this range.

The transmittance spectra calculated by the size-independent EMT and the size-dependent EMT are illustrated in Fig. 4 to show the applicability of the two methods. Both methods predict the existence of the two absorption characteristics successfully. The size-dependent EMT works better than the size-independent EMT at all incident angles. For s polarization, the size-dependent EMT can more accurately predict the width of the absorption peak near 380 nm. However, the transmittance predicted by the two methods is lower than the results calculated by FEM when the wavelength is longer than 500 nm. The two methods have a significant error in predicting the transmittance for p polarization. There is a deviation in the depth of the absorption peak near 600 nm. The deviation keeps increase with the diameter. Moreover, FEM results indicate that the absorption peak undergoes a blueshift with the increase of the incident angle, especially in the nanorods of *d* = 20 nm. Both methods failed to predict the phenomenon. In addition, the two methods are not limited by the Kramers-Kronig relationship. Therefore, the local methods cannot describe the detailed characteristics of nanorod arrays.

Although theoretically speaking, it is more reasonable to use the *k*-dependent effective permittivity to describe the characteristics around 600 nm, the prediction error of the nanorod array with a small diameter is lower. For the nanorod array of *d* = 4 nm, the MSDT for p polarization incident light is reduced to 1.62%, 1.71%, 1.73% and 2.00%. MSDT at all incident angles meets the standard. The results indicate that nanorods of small diameters have the potential to avoid the fluctuations caused by optical nonlocality. The error gradually becomes obvious as the diameter increases.

The extracted effective permittivities are displayed in Fig. 5. The results of the local size-independent EMT, which is used as the original value, are also displayed for comparison. The in-plane permittivity basically follows the oscillator characteristics predicted by EMT. The line widths of the absorption peaks near 400 nm of the three nanorods are slightly different. Therefore, the size of nanorods shows a slight effect on the widths of oscillators. The real part of out-of-plane permittivity decrease from positive to negative, similar to a Drude metal. The position of the absorption peaks gradually shifts toward short wavelengths with the increase of diameter. As a result, the ENZ positions of nanorods of *d* = 4 nm, *d* = 12 nm, *d* = 20 nm are 614 nm, 605 nm, and 596 nm, respectively. The imaginary part of the out-of-plane permittivity decreases rapidly in the short-wavelength region and gradually increases with wavelengths in the long-wavelength region. As proved in previous work, the optical nonlocality mainly works near the ENZ region [38]. When the real part of the out-of-plane permittivity tends to zero, different from the local EMT, the imaginary part appears to fluctuate significantly. The amplitude of the fluctuation increases with the diameter of nanorods. This may indicate that optical nonlocality has a stronger impact on nanorods of larger diameters. As a result, the fitting effect of the transmission spectrum will be worse with the increase of diameter.

## 4. Fabrication of ultrathin nanorods and extraction of effective parameters

The above calculations manifest that the decrease in diameter can reduce the influence of nonlocality. Therefore, we fabricate an ultrathin nanorod array with a diameter of only about 4 nm. The sample is created through lithography-free one-step radio frequency (RF) sputtering deposition. This process is performed via a co-sputtering system (Jsputter8000, ULVAC Inc.). Pure Ag and Al_{2}O_{3} targets are operated together, and an RF is appended for plasma-assisted deposition. The sputter powers of Ag and Al_{2}O_{3} targets are 10 W and 120 W, respectively, with an RF substrate bias of 45 W. The surface free energy of the {111} plane of the silver fcc structure is lower than that of the other planes; therefore, with the assistance of RF, the silver atoms tend to accumulate along the direction perpendicular to the substance. This method has proven to be effective for fabricating ultrathin nanorod arrays of high optical quality [41]. The characteristic of the sample is shown in Fig. 6(a). The transmission spectra for s and p polarization incident light at 20°, 25°, 30°, and 60° are measured by ellipsometer and shown in Fig. 6(b) and Fig. 6(c). The experimental transmittance shows absorption characteristics similar to the calculated results. The positions of the two absorption peaks are about 370 nm and 515 nm, respectively. However, the absorption peak near the position of ENZ is wider than the calculated results, obviously. The difference in transmittance may arise from the limitation of the small diameter to the movement of electrons [30]. The restricted mean free path of electrons will affect the permittivity of silver. The extracting range is set as 270 to 1000 nm to include the two absorption peaks completely.

We get the effective permittivity tensor of the sample through the method above. The transmittance predicted by effective permittivity is shown by the dots in Fig. 6(b) and (c). Although the effective permittivity is *k*-independent, the predicted curve shows great similarity to the measurement results. The MSDTs for s and p polarization incident light are shown in Table 3. Remarkably, the MSDT values for every polarization and all incident angles are lower than 2%. Therefore, it is reasonable to express the transmittance behavior of ultrathin nanorods by the wave-vector independent effective permittivity. The best-matched permittivity is illustrated in Fig. 7. The basic trend of permittivity is consistent with the simulation results. This confirms the reliability of the obtained parameters. The ENZ position is 537 nm, slightly longer than 515 nm where the absorption is strongest. It manifests that the ENZ position can be roughly determined from the position of absorption peak.

## 5. Extraction of *k*-dependent effective permittivity

To obtain more accurate effective permittivity of nanorod arrays of *d* = 12 nm and *d* = 20 nm, we extract a set of effective permittivity that depends upon the incident angle. This model is Kramers-Kronig relationship-satisfied and can obviously reduce the error caused by optical nonlocality. We denote it as *k*-dependent effective permittivity. The above calculation suggests that it is possible to find a set of in-plane permittivity to describe the transmittance spectra for s polarization incident light. Therefore, like previous reports, we suppose the out-of-plane permittivity will vary with the incident angles due to the influence of optical nonlocality [30,50,53]. We keep the extracted in-plane permittivity constant and use the *k*-independent out-of-plane permittivity as the initial value, then extract four sets of out-of-plane permittivities respectively from the transmittance spectra for p polarization incident light at each single incident angle (20°, 25°, 30°, and 60°). Thus, four sets of *k*-dependent effective permittivities are obtained. We find that the transmittance spectra at 20°, 25°, 30° can be expressed by the same set of permittivity. The transmittance spectra at 60° need to be characterized by another set of permittivity. The extracted *k*-dependent out-of-plane permittivity and the predicted transmission spectra are shown in Fig. 8. The MSDT values are shown in Table 4, which are obviously lower than the *k*-independent results at all incident angles. The results indicate that the *k*-dependent effective permittivity is more reasonable for expressing the transmittance spectra of large-diameter nanorod arrays.

The main difference between the two sets of out-of-plane permittivities is that the imaginary part of the permittivity extract from 60° is significantly lower than that obtain from 20°, 25°, and 30° near the ENZ position. This phenomenon is more evident in the nanorod arrays of *d* = 20 nm. This means that optical nonlocality has a greater impact on the transmittance of nanorod arrays of large diameters. The ENZ position also shifts toward short wavelengths with the increasing incident angle. In the nanorod array of *d* = 12 nm, the position shifts from 608 nm to 600 nm. In the nanorod array of *d* = 20 nm, the position shifts from 602 nm to 589 nm. Therefore, the influence of nonlocality is the reduction in the imaginary part and the blueshift of ENZ position in the real part of out-of-plane permittivity during the increasing of incident angle. This may reveal why only the optical properties of small-diameter nanorod arrays can be described by the *k*-independent effective permittivity. In HMMs, the photonic bandgap will undergo a blueshift as the incident angle increases to meet the Bragg condition in general [17,54]. However, the blueshift of the absorption peak of nanorod arrays is stronger than that of the *k*-independent model. A stronger blue shift means a greater phase change as the incident angle changes. Therefore, the main reason why nanorods of large diameter cannot be expressed by effective permittivity may be that the *k*-independent effective permittivity cannot express the phase change of nanorod arrays.

## 6. Conclusion

In this work, we extract the Kramers-Kronig relationship-satisfied in-plane and out-of-plane effective complex permittivity of nanorod arrays with optical nonlocality in the visible range. To quantitative show how to control the effective permittivity by the diameter and incident angle, transmission spectra calculated from 20° to 60° of nanorods with diameters of 4 nm, 12 nm, and 20 nm are selected for calculation. The optical nonlocality leads to a blueshift of the ENZ position in Re(*ɛ*_{⊥}) and a fluctuation in Im(*ɛ*_{⊥}). The blueshift and fluctuation increase with the diameter of nanorods. We demonstrate that transmittance spectra of nanorod array can be expressed by effective permittivity when the diameter is reduced to 4 nm. We also obtain *k*-dependent effective permittivity to provide a more accurate solution for nanorod arrays with large diameters. The *k*-dependent effective permittivity can obviously reduce the error caused by optical nonlocality. The MSDT of nanorods of *d* = 12 nm reduces to less than 2% at all selected incident angles. And the MSDT of nanorods of *d* = 20 nm reduces by more than 1.2%. The results show that when the incident angle increases from 20° to 60°, the ENZ position of Re(*ɛ*_{⊥}) shifts to shorter wavelengths by 8 nm in the nanorod arrays of *d* = 12 nm, and by 13 nm in the nanorod arrays of *d* = 20 nm. This study is significant for the characterization and application of materials with optical nonlocality.

## Funding

National Natural Science Foundation of China (11474078).

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## Data availability

Data underlying the results presented in this paper are available from the corresponding author upon request.

## References

**1. **A. Poddubny, I. Iorsh, P. Belov, and Y. Kivshar, “Hyperbolic metamaterials,” Nat. Photonics **7**(12), 948–957 (2013). [CrossRef]

**2. **G. A. Wurtz, R. Pollard, W. Hendren, G. P. Wiederrecht, D. J. Gosztola, V. A. Podolskiy, and A. V. Zayats, “Designed ultrafast optical nonlinearity in a plasmonic nanorod metamaterial enhanced by nonlocality,” Nat. Nanotechnol. **6**(2), 107–111 (2011). [CrossRef]

**3. **L. H. Nicholls, T. Stefaniuk, M. E. Nasir, F. J. Rodríguez-Fortuo, G. A. Wurtz, and A. V. Zayats, “Designer photonic dynamics by using non-uniform electron temperature distribution for on-demand all-optical switching times,” Nat. Commun. **10**(1), 2967–2968 (2019). [CrossRef]

**4. **A. D. Neira, N. Olivier, M. E. Nasir, W. Dickson, G. A. Wurtz, and A. V. Zayats, “Eliminating material constraints for nonlinearity with plasmonic metamaterials,” Nat. Commun. **6**(1), 7757 (2015). [CrossRef]

**5. **J. Zhou, A. F. Kaplan, L. Chen, and L. J. Guo, “Experiment and theory of the broadband absorption by a tapered hyperbolic metamaterial array,” ACS Photonics **1**(7), 618–624 (2014). [CrossRef]

**6. **P. N. Dyachenko, S. Molesky, A. Y. Petrov, M. StöRmer, T. Krekeler, S. Lang, M. Ritter, Z. Jacob, and M. Eich, “Controlling thermal emission with refractory epsilon-near-zero metamaterials via topological transitions,” Nat. Commun. **7**(1), 11809 (2016). [CrossRef]

**7. **I. A. Kolmychek, A. R. Pomozov, A. P. Leontiev, K. S. Napolskii, and T. V. Murzina, “Magneto-optical effects in hyperbolic metamaterials,” Opt. Lett. **43**(16), 3917–3920 (2018). [CrossRef]

**8. **M. Habib, D. Briukhanova, N. Das, B. C. Yildiz, and H. Caglayan, “Controlling the plasmon resonance via epsilon-near-zero multilayer metamaterials,” Nanophotonics **9**(11), 3637–3644 (2020). [CrossRef]

**9. **R. Yu, R. Alaee, R. W. Boyd, and F. J. Abajo, “Ultrafast topological engineering in metamaterials,” Phys. Rev. Lett. **125**(3), 037403 (2020). [CrossRef]

**10. **F. Wu, “Achieving highly efficient and broad-angle polarization beam filtering using epsilon-near-zero metamaterials mimicked by metal-dielectric multilayers,” J. Appl. Phys. **123**(9), 093103 (2018). [CrossRef]

**11. **C. H. Xue, Y. Ding, H. T. Jiang, Y. Li, Z. Wang, Y. Zhang, and H. Chen, “Dispersionless gaps and cavity modes in photonic crystals containing hyperbolic metamaterials,” Phys. Rev. B **93**(12), 125310 (2016). [CrossRef]

**12. **T. Li and J. B. Khurgin, “Hyperbolic metamaterials: Beyond the effective medium theory,” Optica **3**(12), 1388 (2016). [CrossRef]

**13. **J. Sukham, O. Takayama, M. Mahmoodi, S. Sychev, A. Bogdanov, S. H. Tavassoli, A. V. Lavrinenkoa, and R. Malureanu, “Investigation of effective media applicability for ultrathin multilayer structures,” Nanoscale **11**(26), 12582–12588 (2019). [CrossRef]

**14. **V. Popov and A. Novitsky, “Advanced effective medium approximation for subwavelength multilayers to overcome Maxwell-Garnet approach breakdown,” in 10th International Congress on Advanced Electromagnetic Materials in Microwaves and Optics, 2016, pp. 289–291.

**15. **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**(8), 085428 (2016). [CrossRef]

**16. **V. Popov, A. V. Lavrinenko, and A. Novitsky, “Surface waves on multilayer hyperbolic metamaterials: Operator approach to effective medium approximation,” Phys. Rev. B **97**(12), 125428 (2018). [CrossRef]

**17. **W. Feng, G. Lu, C. Xue, H. Jiang, Z. Guo, M. Zheng, C. Chen, G. Du, and C. Hong, “Experimental demonstration of angle-independent gaps in one-dimensional photonic crystals containing layered hyperbolic metamaterials and dielectrics at visible wavelengths,” Appl. Phys. Lett. **112**(4), 041902 (2018). [CrossRef]

**18. **M. A. Noginov, Y. A. Barnakov, G. Zhu, T. Tumkur, H. Li, and E. E. Narimanov, “Bulk metamaterial with hyperbolic dispersion,” Appl. Phys. Lett. **94**(15), 151105 (2009). [CrossRef]

**19. **T. Tumkur, Y. Barnakov, S. T. Kee, M. A. Noginov, and V. Liberman, “Permittivity evaluation of multilayered hyperbolic metamaterials: Ellipsometry vs. reflectometry,” J. Appl. Phys. **117**(10), 103104 (2015). [CrossRef]

**20. **S. Hao, D. Lu, B. Vansaders, J. J. Kan, and Z. Liu, “Anomalously weak scattering in metal-semiconductor multilayer hyperbolic metamaterials,” Phys. Rev. X **5**(2), 021021 (2015). [CrossRef]

**21. **O. Y. Yermakov, D. V. Permyakov, F. V. Porubaev, P. A. Dmitriev, A. K. Samusev, I. V. Iorsh, R. Malureanu, A. V. Lavrinenko, and A. A. Bogdanov, “Effective surface conductivity of optical hyperbolic metasurfaces: from far-field characterization to surface wave analysis,” Sci. Rep. **8**(1), 14135 (2018). [CrossRef]

**22. **C. Zhang, N. Hong, C. Ji, W. Zhu, X. Chen, A. Agrawal, Z. Zhang, T. E. Tiwald, S. Schoeche, J. N. Hilfiker, L. J. Guo, and H. J. Lezec, “Robust extraction of hyperbolic metamaterial permittivity using total internal reflection ellipsometry,” ACS Photonics **5**(6), 2234–2242 (2018). [CrossRef]

**23. **J. Dilts, C. Hong, A. Siahmakoun, M. Syed, and H. Alisafaee, “Low-MSE extraction of permittivity in optical hyperbolic metamaterials,” Opt. Lett. **44**(17), 4303–4306 (2019). [CrossRef]

**24. **B. Gompf, J. Braun, T. Weiss, H. Giessen, M. Dressel, and U. Hübner, “Periodic nanostructures: Spatial dispersion mimics chirality,” Phys. Rev. Lett. **106**(18), 185501 (2011). [CrossRef]

**25. **M. A. Gorlach, S. B. Glybovski, A. A. Hurshkainen, and P. A. Belov, “Giant spatial-dispersion-induced birefringence in metamaterials,” Phys. Rev. B **93**(20), 201115 (2016). [CrossRef]

**26. **M. A. Gorlach and P. A. Belov, “Effect of spatial dispersion on the topological transition in metamaterials,” Phys. Rev. B **90**(11), 115136 (2014). [CrossRef]

**27. **M. A. Shapiro, G. Shvets, J. R. Sirigiri, and R. J. Temkin, “Spatial dispersion in metamaterials with negative dielectric permittivity and its effect on surface waves,” Opt. Lett. **31**(13), 2051–2053 (2006). [CrossRef]

**28. **J. M. Mcmahon, S. K. Gray, and G. C. Schatz, “Nonlocal optical response of metal nanostructures with arbitrary shape,” Phys. Rev. Lett. **103**(9), 097403 (2009). [CrossRef]

**29. **R.-L. Chern and D. Han, “Nonlocal optical properties in periodic lattice of graphene layers,” Opt. Express **22**(4), 4817–4829 (2014). [CrossRef]

**30. **B. Janaszek, M. Kieliszczyk, A. Tyszka-Zawadzka, and P. Szczepański, “Influence of nonlocality on transmittance and reflectance of hyperbolic metamaterials,” Crystals **10**(7), 577 (2020). [CrossRef]

**31. **B. Janaszek and P. Szczepański, “Effect of nonlocality in spatially uniform anisotropic metamaterials,” Opt. Express **28**(10), 15447–15458 (2020). [CrossRef]

**32. **S. Gong, M. Hu, Z. Wu, H. Pan, H. Wang, K. Zhang, R. Zhong, J. Zhou, T. Zhao, D. Liu, W. Wang, C. Zhang, and S. Liu, “Direction controllable inverse transition radiation from the spatial dispersion in a graphene-dielectric stack,” Photonics Res. **7**(10), 1154–1160 (2019). [CrossRef]

**33. **J. Benedicto, R. Pollès, C. Ciracì, E. Centeno, D. R. Smith, and A. Moreau, “Numerical tool to take nonlocal effects into account in metallo-dielectric multilayers,” J. Opt. Soc. Am. A **32**(8), 1581–1588 (2015). [CrossRef]

**34. **L. Sun, J. Gao, and X. Yang, “Giant optical nonlocality near the Dirac point in metal-dielectric multilayer metamaterials,” Opt. Express **21**(18), 21542–21555 (2013). [CrossRef]

**35. **K. J. Lee, Y. U. Lee, F. Fages, J. C. Ribierre, J. W. Wu, and A. D’Aléo, “Blue-shifting intramolecular charge transfer emission by nonlocal effect of hyperbolic metamaterials,” Nano Lett. **18**(2), 1476–1482 (2018). [CrossRef]

**36. **B. Wells, Z. A. Kudyshev, N. Litchinitser, and V. A. Podolskiy, “Nonlocal effects in transition hyperbolic metamaterials,” ACS Photonics **4**(10), 2470–2478 (2017). [CrossRef]

**37. **M. G. Silveirinha, “Nonlocal Homogenization Model for a Periodic Array of Epsilon-Negative Rods,” Phys. Rev. E **73**(4), 046612 (2006). [CrossRef]

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

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

**40. **L. V. Alekseyev, E. E. Narimanov, T. Tumkur, H. Li, Y. A. Barnakov, and M. A. Noginov, “Uniaxial epsilon-near-zero metamaterial for angular filtering and polarization control,” Appl. Phys. Lett. **97**(13), 131107 (2010). [CrossRef]

**41. **J. Gao, X. Wu, Q. Li, S. Du, F. Huang, L. Liang, H. Zhang, F. Zhuge, H. Cao, and Y. Song, “Template-free growth of well-ordered silver nano forest/ceramic metamaterial films with tunable optical responses,” Adv. Mater. **29**(16), 1605324 (2017). [CrossRef]

**42. **S. Peruch, A. Neira, G. A. Wurtz, B. Wells, V. A. Podolskiy, and A. V. Zayats, “Geometry defines ultrafast hot-carrier dynamics and Kerr nonlinearity in plasmonic metamaterial waveguides and cavities,” Adv. Opt. Mater. **5**(15), 1700299 (2017). [CrossRef]

**43. **H.-J. Hagemann, W. Gudat, and C. Kunz, “Optical constants from the far infrared to the x-ray region: Mg, Al, Cu, Ag, Au, Bi, C, and A1_{2}O_{3},” J. Opt. Soc. Am. **65**(6), 742–744 (1975). [CrossRef]

**44. **A. P. Bradford, G. Hass, and M. McFarland, “Optical properties of evaporated magnesium oxide films in the 0.22-8micro wavelength region,” Appl. Opt. **11**(10), 2242–2244 (1972). [CrossRef]

**45. **R. Atkinson, W. R. Hendren, G. A. Wurtz, W. Dickson, A. V. Zayats, P. Evans, and R. J. Pollard, “Anisotropic optical properties of arrays of gold nanorods embedded in alumina,” Phys. Rev. B **73**(23), 235402 (2006). [CrossRef]

**46. **B. Gompf, B. Krausz, B. Frank, and M. Dressel, “*k*-dependent optics of nanostructures: Spatial dispersion of metallic nanorings and split-ring resonators,” Phys. Rev. B **86**(7), 075462 (2012). [CrossRef]

**47. **V. M. Agranovich and V. Ginzburg, * Crystal optics with spatial dispersion and excitons*, (Springer, 1984), Chap. 1.

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

**49. **H. Y. Xie, M. Y. Ng, and Y. C. Chang, “Analytical solutions to light scattering by plasmonic nanoparticles with nearly spherical shape and nonlocal effect,” J. Opt. Soc. Am. A **27**(11), 2411–2422 (2010). [CrossRef]

**50. **T. Geng, S. Zhuang, J. Gao, and X. Yang, “Nonlocal effective medium approximation for metallic nanorod metamaterials,” Phys. Rev. B **91**(24), 245128 (2015). [CrossRef]

**51. **B. Johs and J. S. Hale, “Dielectric function representation by B-splines,” Phys. Status Solidi A **205**(4), 715–719 (2008). [CrossRef]

**52. **K. Madsen, H.B. Nielsen, and O. Tingleff, Methods for Non-Linear Least Squares Problems, 2nd ed.

**53. **S. I. Pekar, “The theory of electromagnetic waves in a crystal in which excitons are produced,” J. Phys. Chem. Solids **5**(1-2), 11–22 (1958). [CrossRef]

**54. **J. N. Winn, Y. Fink, S. Fan, and J. D. Joannopoulos, “Omnidirectional reflection from a one-dimensional photonic crystal,” Opt. Lett. **23**(20), 1573–1575 (1998). [CrossRef]