Enhanced nonlocal effects in metamaterials with moderate-index inclusions

TL;DR

This paper explores multilayered metamaterials with moderate-index inclusions, revealing significant nonlocal effects and novel dispersion phenomena that could advance optical material engineering.

Contribution

It demonstrates non-resonant enhancement of spatial dispersion effects in moderate-index inclusions, enabling new dispersion regimes and topological phase transitions.

Findings

Enhanced nonlocal effects comparable to local response

Identification of a novel configuration with disconnected isofrequency contours

Potential for practical optical applications without high-index inclusions

Abstract

We investigate a class of multilayered metamaterials characterized by moderate-index inclusions and low average permittivity. Via first-principle calculations, we show that in such scenario first- and second-order spatial dispersion effects may exhibit a dramatic and non-resonant enhancement, and may become comparable to the local response. Their interplay gives access to a wealth of dispersion regimes encompassing additional extraordinary waves and topological phase transitions. In particular, we identify a novel configuration featuring bound and disconnected isofrequency contours. Since they do not rely on high-index inclusions, our proposed metamaterials may constitute an attractive and technologically viable platform for engineering nonlocal effects in the optical range.

Equations

$$|{\epsilon }_B|$$

$$|\overline{\epsilon (x)}|$$

$${\overline{D}}_x$$

$${\overline{D}}_z$$

$${\overline{B}}_y$$

$${ϵ}_{||}$$

$$\kappa$$

$$\delta$$

$$\delta q_z^4-\left({ϵ}_{\perp }+{\kappa }^2+\delta {ϵ}_{||}\right)q_z^2-{ϵ}_{||}{q}_x^2+{ϵ}_{||}{ϵ}_{\perp }=0.$$

$${\partial }_z{E}_x-{\partial }_x{E}_z$$

$${\partial }_z{H}_y$$

$${\partial }_x{H}_y$$

$$D_x$$

$$B_y$$

$$ϵ$$

$${ϵ}^{-1}$$

$$E_x$$

$$E_z$$

$$H_y$$

$${\partial }_z{\overline{E}}_x-{\partial }_x{\overline{E}}_z=i{\mu }_0\omega {\overline{H}}_y,$$

$$\overline{{ϵ}^{-1}}{\partial }_z{\overline{H}}_y+\overline{\widetilde{{ϵ}^{-1}}{\partial }_z{\tilde{H}}_y}=i\omega {ϵ}_0{\overline{E}}_x,$$

$${\partial }_x{\overline{H}}_y=-i\omega {ϵ}_0\left(\overline{ϵ}{\overline{E}}_z+\overline{\widetilde{ϵ}{\tilde{E}}_z}\right),$$

$${\partial }_z{\tilde{E}}_x-{\partial }_x{\tilde{E}}_z-i{\mu }_0\omega {\tilde{H}}_y=0,$$

$$i\omega {ϵ}_0{\tilde{E}}_x-{ϵ}^{-1}{\partial }_z{\tilde{H}}_y+\overline{\widetilde{{ϵ}^{-1}}{\partial }_z{\tilde{H}}_y}=\widetilde{{ϵ}^{-1}}{\partial }_z{\overline{H}}_y,$$

$${\partial }_X{\tilde{H}}_y=-i\omega {ϵ}_0\widetilde{ϵ}{\overline{E}}_z.$$

$$\begin{gathered} \widetilde{ϵ}=\sum _{n=-\infty \\ n\ne 0}^{+\infty }{A}_n{e}^{in{k}_{0}X},\widetilde{{ϵ}^{-1}}=\sum _{n=-\infty \\ n\ne 0}^{+\infty }{B}_n{e}^{in{k}_{0}X}, \end{gathered}$$

$${\partial }_z{\overline{E}}_x-{\partial }_x{\overline{E}}_z$$

$${\partial }_z{\overline{H}}_y$$

$${\partial }_x{\overline{H}}_y$$

$${\overline{D}}_x$$

$${\overline{D}}_z$$

$${\overline{B}}_y$$

$${ϵ}_{||}$$

$$\kappa$$

$$\beta$$

$${\overline{D}}_x$$

$${\overline{D}}_z$$

$${\overline{B}}_y$$

$$\delta =\frac{\beta -{\kappa }^2}{{ϵ}_{||}}.$$

$${\overline{D}}_x$$

$${\overline{D}}_z$$

$${\overline{B}}_y$$

$${ϵ}_{||}^S$$

Full text

Metamaterials are artificial composites of dielectric and/or metallic inclusions in a host medium, which can be engineered so as to exhibit unconventional and/or tailored effective electromagnetic responses. As such, they constitute a unique platform for attaining novel optical mechanisms such as negative refraction, super- and hyper-lensing, invisibility cloaking, etc. Cai . In the long-wavelength limit, the metamaterial optical response is typically described in terms of macroscopic phenomenological parameters, such as effective permittivity and permeability Smith . However, if the inclusions are metallic and/or their size is not electrically small, nonlocal effects (i.e., spatial dispersion) may need to be accounted for, e.g., via the appearance of spatial derivatives of the fields in the effective constitutive relationships, or via wavevector-dependent constitutive parameters Landau . The reader is referred to Silveirinha ; Elser ; Alu ; Chebykin1 ; Chebykin2 for a representative sampling of nonlocal homogenization approaches available in the topical literature. It is worth stressing that two cornerstones in metamaterial science, namely, artificial electromagnetic chirality and optical magnetism, are manifestations of first- and second-order spatial dispersion, respectively Serd ; Ciatt_1 . Artificial electromagnetic chirality (due to 3-D, 2-D and 1-D geometrical chirality Rizza_1 ) yields interesting phenomena, such as giant optical activity, asymmetric transmission, and negative refractive index Li . Artificial or optical magnetism may give rise to negative refractive index and backward-wave propagation Smith_2000 .

It is worth noting that second-order spatial dispersion is not fully equivalent to artificial magnetism Gorlach , and this aspect has not been exploited to its fullest extent in applied science. Even though second-order and non-magnetic nonlocal effects are generally considered detrimental for several applications Taming , they have been shown to support intriguing phenomena such as propagation of additional extraordinary waves Orlov1 ; Zayats and topological transitions Gorlach ; Orlov2 ; Mirmo . Within this context, harnessing these nonlocal effects currently represents one of the grand challenges of metamaterial science.

Spatial dispersion is essentially ruled by the electrical size $\eta=\Lambda/\lambda$, with $\Lambda$ denoting a characteristic metamaterial scale (unit-cell size) and $\lambda$ the vacuum wavelength. Therefore, unlike conventional materials, metamaterials offer a natural way for enhancing nonlocality by tailoring the spatial distribution of the inclusions. There are several strategies to modify the hierarchy of spatial-dispersion orders: i) by optimizing shape and size of the inclusions Saenz , ii) by exploiting plasmonic resonances Leon , iii) by resorting to extreme regimes where permittivity and/or permeability are close to zero Rizza_1 , and iv) by embedding high-index dielectric inclusions, i.e., with relative permittivity $|\varepsilon_{I}|\sim 1/\eta^{2}$ Felbacq ; Merlin . The latter strategy has led to the development of low-loss, all-dielectric metamaterials with electric and magnetic responses Krasnok .

Here, we propose a strategy to enhance both first- and second-order nonlocal effects, so that their interplay may become comparable to the local response, without resorting to plasmonic resonances. Such interplay is remarkably unusual in non-resonant scenarios, since first- and second-order spatial-dispersion effects typically differ by orders of magnitude. Moreover, their symmetry properties are also different, since first-order nonlocality is ruled by geometric chirality, whereas second-order nonlocality is not.

We show that the above strategy can be implemented in multilayered metamaterials characterized by moderate-index inclusions (MII), i.e., with relative permittivity $|\varepsilon_{I}|\sim 1/\eta$, in conjunction with a relatively low average relative permittivity ($\ll 1/\eta$). Via first-principle derivations, we show that transverse magnetic (TM) waves are ruled by a biquadratic, highly tunable dispersion relationship that yields a broad variety of phenomena, encompassing additional extraordinary waves and topological phase transitions. Besides the well-known hyperbolic, elliptic and mixed regimes, we identify a topologically nontrivial dispersion regime characterized by bounded and multiply-connected isofrequency contours accompanied by additional extraordinary waves. The proposed metamaterials, which combine a fairly simple geometrical structure with relaxed (moderate-index) requirements on the material constituents, may provide a versatile and technologically viable platform for dispersion engineering, especially attractive in spectral regions (e.g., optical range) where extremely high-index materials are not available.

With reference to the schematic in Fig. 1(a), we consider a 1-D multilayered metamaterial, periodically stratified along the $x$-axis, with a unit cell of thickness $\Lambda$ consisting of two (possibly composite) regions that we generically indicate as “background” and “inclusion”. All fields and quantities are assumed as independent of $y$, and the metamaterial is described by a piecewise continuous relative permittivity distribution $\varepsilon(x)$, denoted by $\varepsilon_{B}$ and $\varepsilon_{I}$ within the background and the inclusion regions, respectively. In the presence of a monochromatic electromagnetic field (with suppressed $e^{-i\omega t}$ time-dependence) whose vacuum wavelength $\lambda$ is much larger than the unit-cell period, the small ratio $\eta=\Lambda/\lambda$ rules the homogenization of the metamaterial response. As anticipated, we are interested in the MII regime characterized by the basic requirements:

[TABLE]

with the overline indicating the average over the unit cell. As detailed in Supplement 1, a rigorous, multi-scale analysis of TM-wave ($y$-directed magnetic field) propagation, in the long-wavelength limit $\eta\ll 1$ and the MII regime, shows that the leading-order metamaterial response is governed by the effective constitutive relationships

[TABLE]

with $\epsilon_{0}$, $\mu_{0}$ and $k_{0}=\omega\sqrt{\epsilon_{0}\mu_{0}}=2\pi/\lambda$ denoting the vacuum permittivity, permeability and wavenumber, respectively, and

[TABLE]

where the summations are intended over all integers $n,m\neq 0$, with $n\neq-m$, and $A_{n}/\eta$ and $B_{n}$ are the Fourier coefficients of the functions $\varepsilon(x)$ and $\varepsilon^{-1}(x)$, respectively (with $|A_{n}|,\>|B_{n}|\ll 1/\eta$). Equations (3) and (33) summarize the main result of this study, and they show that the MII regime, in addition to the local response, exhibits both first- and second-order spatial dispersion. By using the Serdyukov-Fedorov transformation method Serd , it can be shown that the obtained first-order spatial dispersion is equivalent to 1-D electromagnetic chirality Rizza_1 , and that the second-order spatial dispersion is non-magnetic (i.e., there is no transformation yielding a non-vanishing magnetic permeability). It is crucial to stress that Eqs. (3) do not contain the small parameter $\eta$. From the physical viewpoint, this yields an unprecedented regime where local response and nonlocal effects may interplay at comparable levels. By contrast, in media with shallow dielectric modulation, first- and second-order spatial-dispersion contributions scale as $\sim\eta$ and $\sim\eta^{2}$, respectively, and hence are too weak to appear in the leading-order term of the effective response Merlin .

Our proposed MII regime in Eqs. (1) forbids the onset of resonances that would produce a dominant local response, and, at the same time, it triggers an anomalous enhancement of nonlocal effects which strengthens (in different ways) the relative weight of standard first- and second-order spatial-dispersion contributions, while leaving higher-order spatial-dispersion contributions at a negligibly weak level. To gain some physical insight into the underlying mechanisms, we start analyzing the local response. Due to the requirements in Eqs. (1), the average of $\varepsilon^{-1}(x)$ is dominated by the background contribution $\varepsilon_{B}$, so that $|\epsilon_{||}|\ll 1/\eta$ [see Eq. (6)]. This is a key aspect characterizing the MII regime since, in the presence of a uniformly large dielectric spatial modulation, $|\epsilon_{||}|$ can be exceedingly large so as to yield superlensing effects Rizza_2 which are unaffected by spatial dispersion. The requirements in Eqs. (1), attainable by suitably tailoring the inclusion dielectric profile (see our discussion and examples below), imply that the first contribution $\overline{\varepsilon(x)}$ to $\epsilon_{\perp}$ in Eq. (6) is not large. Moreover, the Fourier coefficients of $\varepsilon(x)$ (with $n\neq 0$) are $A_{n}/\eta$ with $|A_{n}|\ll 1/\eta$, and hence $|\epsilon_{\perp}|\ll 1/\eta$. In other words, even if the dielectric modulation is large, the MII regime is characterized by local effective permittivities that are sufficiently small so as not to overshadow the spatial-dispersion effects. Considering now the nonlocal response, the requirements in Eqs. (1) imply that all Fourier coefficients of $\varepsilon^{-1}(x)$ are dominated by the background permittivity, and hence $|B_{n}|\ll 1/\eta$. As a consequence, from Eqs. (7) and (8), the strength of the nonlocal response may become comparable to the local one (see also Supplement 1 for a heuristic explanation of the enhancement effects). It is also worth emphasizing that, although the MII requirements in Eqs. (1) do not explicitly provide indications on the permittivity signs, the most significant spatial-dispersion enhancements [i.e., respectable values of $\kappa$ and $\delta$ in Eqs. (33)] are obtained in configurations mixing positive- and negative-permittivity constituents. From the physical viewpoint, this is not surprising, since these configurations may support surface-plasmon-polaritons propagating at the layer interfaces, whose coupling is known to yield strong nonlocal effects Kidwai .

In order to illustrate the physical implications of the strong MII nonlocal response, we analyze the dispersion relationship of TM waves $H_{y}=U_{y}e^{ik_{0}(q_{x}x+q_{z}z)}$, which readily follows from the constitutive relationships in Eqs. (3), viz.,

[TABLE]

A first important consequence of second-order spatial-dispersion (represented by $\delta$) is that such equation is biquadratic. Moreover, the four coefficients in Eq. (9) depend on the four independent effective parameters in Eqs. (33) which, in turn, may assume comparable values. This implies that, in principle, every possible dispersion regime arising from Eq. (9) is accessible.

Figures 1(b)–1(f) qualitatively illustrate, in terms of isofrequency contours, the five possible and topologically different dispersion regimes that can be attained: hyperbolic-like I and II, mixed, elliptic-like I and II, respectively. It is remarkable that the relatively simple class of MII metamaterials is able to support such a variety of dispersion regimes, comprising both bounded [Figs. 1(e) and 1(f)] and unbounded [Figs. 1(b), 1(c) and 1(d)] propagating wave spectra whose isofrequency contours can be either connected [Fig. 1(e)] or disconnected [Figs. 1(b), 1(c), 1(d) and 1(f)]. Elliptic-hyperbolic mixed regimes such as in Fig. 1(d) have been predicted and experimentally observed Zayats ; Gorlach ; Mirmo , together with additional-wave manifestations, in different metamaterial configurations where second-order spatial dispersion plays a key role in the absence of electromagnetic chirality. To the best of our knowledge, the elliptic-like II regime in Fig. 1(f) has never been predicted in 1-D metamaterials; its distinctive feature is the presence of two propagating waves with the same transverse wavenumber (i.e., an additional extraordinary wave) within a limited $q_{x}$ spectral region.

To validate and calibrate the MII effective model in Eqs. (3), we compare its predictions with rigorous Bloch theory Joannopoulos and full-wave numerical simulations Comsol (see Supplement 1 for details). As a first illustrative example, we consider a unit cell comprising three homogeneous layers, two of which (labeled as $1$ and $2$) constituting the inclusion, and the third one (labeled as $3$) representing the background [see Fig. 2(a)]. We consider a parameter configuration with $\eta=0.1$, layer relative permittivities $\epsilon_{1}=3.8$, $\epsilon_{2}=-4.2$ and $\epsilon_{3}=-0.2$, and filling fractions $f_{1}=f_{2}=0.4$ and $f_{3}=0.2$. Such a design fulfills the MII requirements in Eqs. (1), since $|\varepsilon_{B}|=0.2\ll 1/\eta$, $|\varepsilon_{I}|\simeq 4\sim 1/\eta$ and $|\overline{\varepsilon(x)}|=|f_{1}\epsilon_{1}+f_{2}\epsilon_{2}+f_{3}\epsilon_{3}|=0.2\ll 1/\eta$. For the chosen parameter configuration, we obtain from Eqs. (33) $\epsilon_{\perp}=-9.2\cdot 10^{-2}$, $\epsilon_{||}=-1$, $\kappa=-0.41$ and $\delta=3.4\cdot 10^{-3}$. This yields the elliptic-like II dispersion regime in Fig. 1(f), where additional waves are expected, as sketched in Fig. 2(a). Figure 2(b) compares the isofrequency contours computed from Eqs. (9) with the Bloch-theory Joannopoulos solution. A fairly good agreement is evident, as well as the occurrence of an additional extraordinary wave, which is accurately captured by the MII effective model. Note also that $\overline{\varepsilon(x)}=-0.2$ and $\left[\overline{\varepsilon^{-1}(x)}\right]^{-1}=-1$, so that the standard (local) effective-medium theory (EMT) Cai predicts that plane-wave propagation is forbidden (evanescent) in the medium, in stark contrast with the above evidence. It is worth stressing that, in contrast to high-index-inclusion metamaterials Felbacq , such regime dominated by spatial dispersion is accessed here by means of moderate values of the layer permittivities. As a further validation check, we also analyzed via full-wave simulations a MII metamaterial slab of thickness $L=10$ $\mu$m, with the above parameters and period $\Lambda=0.1$ $\mu$m, under normally incident (along $z$) plane-wave illumination with wavelength $\lambda=1$ $\mu$m [see Fig. 2(a)]. Figure 2(c) shows the magnitude of the average field distribution $|\overline{H}_{y}(z)|$, whose standing-wave pattern reveals the occurrence of different spatial harmonics. This is more evident in the corresponding spatial spectrum (magnitude) shown in Fig. 2(d), which displays four distinct peaks corresponding to two forward and two backward propagating waves, thereby proving the actual excitability of additional extraordinary waves in a finite-thickness slab. For a more quantitative assessment, the normalized wavenumbers $q_{z}$ corresponding to the spectral peaks are superimposed (circle markers) on the isofrequency contours in Fig. 2(b), showing a good agreement with those predicted by the above description of an unbounded MII metamaterial.

The existence of various different dispersion regimes implies that MII metamaterials may also host unusual topological transitions. As a second representative example, we consider a three-layer unit-cell of period $\Lambda=0.05$ $\mu$m, with material constituents chosen as silver (labeled as 1, and playing the role of the inclusion), silicon dioxide (SiO2) and air (labeled as 2,3 respectively, and playing the role of the background), with filling fractions $f_{1}=0.15$, $f_{2}=0.33$, and $f_{3}=0.52$. For silver, we utilize the Drude model $\epsilon_{1}=\epsilon_{b}-\omega_{p}^{2}/(\omega^{2}+i\alpha\omega)$ (with $\epsilon_{b}=5.26$, $\alpha=7.06\cdot 10^{13}$ s*-1*, $\omega_{p}=1.45\cdot 10^{16}$ s*-1* JC ), whereas for SiO2 ($\epsilon_{2}$) we employ an empirical Sellmeier equation Malitson . Figure 3 shows the complex isofrequency contours (real and imaginary parts of $q_{z}$ as a function of $q_{x}$), at three representative wavelengths $\lambda=0.5,0.48,0.46\>\mu m$ for which the MII requirements in Eqs. (1) are fulfilled. Once again, the agreement between the MII effective model and Bloch theory is remarkably good at any wavelength, whereas the EMT provides an inadequate description, as it captures only a portion around $q_{x}=0$ of a single connected component of the isofrequency contours. Note that at, $\lambda=0.5\>\mu m$, the dispersion is of hyperbolic-like I type, since the additional branches [lower in Fig. 3(a), and upper in Fig. 3(b)] are generated by the presence of loss (they disappear in the lossless limit). On the other hand, at $\lambda=0.46\>\mu m$ [Figs. 3(e) and 3(f)], a mixed regime is observed, and this indicates that a topological transition between the two regimes occurs. As also observed in Ref. Orlov2 , such transition is accompanied by an intermediate state at $\lambda=0.48\>\mu m$ [see Figs. 3(c) and 3(d)] where the two dispersion branches of the mixed regime merge at $q_{x}=0$, thereby yielding a wavevector degeneracy and an associated Dirac point.

In conclusion, we have introduced a class of 1-D MII metamaterials that may exhibit an anomalous electrodynamic behavior driven by nonlocal effects. By comparison with other strategies for spatial-dispersion enhancement, our proposed framework relaxes the requirements on the refractive-index of the inclusions, by affording moderate values. In essence, the specific structure of the MII requirements in Eqs. (1) crucially alters the hierarchy of spatial-dispersion orders, leading to a previously unexplored interplay between first- and second-order nonlocal effects, whose strength may become comparable to the local response. This in turns leads to a variety of different dispersion regimes supporting exotic phenomena such as topological phase transitions and additional extraordinary waves. In view of their technological viability and potential versatility, we believe that MII multilayered metamaterials, and their possible higher-dimensional extensions, may play a pivotal role for conceiving novel platforms for dispersion-engineering applications over a broad electromagnetic spectrum, including the optical range.

Supplementary Material

.1 Multiscale Homogenization Theory

Referring to the schematic in Fig. 1(a), we consider a 1-D multilayered metamaterial with a unit-cell of spatial period $\Lambda$, described by a periodic relative permittivity distribution $\varepsilon(x)=\varepsilon(x+\Lambda)$. Moreover, we assume monochromatic, transverse-magnetic (TM) plane-wave propagation, with relevant field components $E_{x}(x,z)$, $E_{z}(x,z)$, $H_{y}(x,z)$ satisfying the source-free Maxwell’s equations

[TABLE]

where

[TABLE]

and all symbols are already defined in the main text. To investigate the long-wavelength limit $\Lambda\ll\lambda$, and to obtain an effective medium description, we introduce the dimensionless small parameter $\eta=\Lambda/\lambda\ll 1$ and the fast-scale coordinate $X=x/\eta$, and we subsequently treat this coordinate as independent from the slow-scale variables $x$ and $z$. As anticipated in the main text, we consider a moderate-index inclusion (MII) parameter configuration fulfilling the requirements in Eqs. (1). In particular, the requirements in Eqs. (1a) imply that $|\varepsilon^{-1}(x)|\ll 1/\eta$ and $\varepsilon-\overline{\varepsilon(x)}\sim 1/\eta$. Accordingly, we address the homogenization of a specific profile $\epsilon(X)=\varepsilon(x)$ (depending only on the fast-scale coordinate $X$) and characterized by the properties

[TABLE]

where the tilde indicates the rapidly varying zero-mean residual. In order to prevent the occurrence of secular terms in the limit $\eta\rightarrow 0$, we consider

[TABLE]

where $\overline{A}$ and $\tilde{A}$ are the slowly and rapidly varying contributions to the field component $A$, respectively. This constitutes the basic Ansatz of the two-scale homogenization method detailed hereafter, where each component $A$ depends on the fast-scale and slow-scale variables, and each electromagnetic field component in Eqs. (17) is a suitable Taylor expansion up to the first order in the parameter $\eta$. By substituting Eqs. (17) and (15) into Eqs. (10), and separating the fast and slow contributions, we obtain

[TABLE]

and

[TABLE]

By introducing the Fourier series

[TABLE]

we solve Eqs. (23), and thereby obtain the rapidly varying contributions in terms of the slow (average) components. By substituting the expressions of the rapidly varying contributions $\widetilde{E}_{x}$, $\widetilde{E}_{z}$ and $\widetilde{H}_{y}$ [derived by Eqs. (23)] into Eqs. (20), we then obtain the equations describing the dynamics of the average fields in the long-wavelength limit. These equations can be rearranged in the form of Maxwell’s equations

[TABLE]

where the effective constitutive relationships are given by

[TABLE]

In Eqs. (30), $Z_{0}=\sqrt{\mu_{0}/\epsilon_{0}}$ denotes the vacuum characteristic impedance, and we have introduced the effective constitutive parameters

[TABLE]

By using Eq. (28), the constitutive relations in Eqs. (30) can be recast in the standard Landau-Lifshitz form Landau , viz.,

[TABLE]

where

[TABLE]

Equations (36) are those reported as Eqs. (2) in the main text. As shown by Eqs. (36), the electromagnetic response is affected by first- and the second-order spatial dispersion, i.e., the first- and second-order derivatives of the electric field describe the medium nonlocality. By using the Serdyukov-Fedorov transformation method Serd , it can be shown that the first-order spatial-dispersion contribution is equivalent to a reciprocal magneto-electric coupling (extrinsic chiral response) whereas the second-order spatial dispersion contribution cannot be cast as a magnetic term. It is important to emphasize that first- and second-order spatial dispersion terms in Eqs. (36) are not weighted by the small parameter $\eta$, which indicates an anomalous enhancement.

.2 Heuristic Interpretation of the Anomalous Spatial-Dispersion Enhancement

The anomalous enhancement of first- and second-order spatial-dispersion effects can be heuristically grasped from the effective medium response of a multilayered metamaterial with shallow dielectric modulation. Following the approach reported in Ref. Ciatt_1 , we obtain the effective constitutive relationships for TM waves propagating in a multilayered metamaterial, viz.,

[TABLE]

The constitutive relationships in the above scenario resemble those in Eqs. (2) [and Eqs. (36)], but without some second-order spatial-dispersion contributions (i.e., $\delta_{a}=0$ and $\delta_{b}=0$), and with parameters (identified by the superscript S)

[TABLE]

where $a_{n}$, $b_{n}$ are the Fourier coefficients of $\varepsilon(x)$ and $\varepsilon^{-1}(x)$, respectively. For this configuration, spatial dispersion is weak since $|\kappa^{S}|\sim\eta\ll 1$, and $|\delta^{S}|\sim|\delta^{S}_{a}|\sim|\delta^{S}_{b}|\sim\eta^{2}\ll 1$ in view of the low dielectric contrast. The crucial observation to make is that if, in Eq. (44) we formally replace $a_{n}$ with $A_{n}/\eta$ and $b_{n}$ with $B_{n}$, we fully recover the expression of $\kappa$ in Eq. (3a). This provides an insightful interpretation of the MII mechanism, suggesting that it effectively elevates the first-order spatial-dispersion effects from the $\eta$ order to the leading one. This observation also sheds light on the enhancement of second-order spatial-dispersion which, in a low-dielectric-contrast medium, is very weak ($\sim\eta^{2}$) with contributions containing either two $a_{n}$ or two $b_{n}$ terms. The above illustrated mechanism of elevating to the leading order the terms $a_{n}=A_{n}/\eta$ (for $n\neq 0$) likewise elevates the contributions containing two $a_{n}$ terms, and this explains their appearance in Eqs. (3). On the other hand, the terms proportional to $\delta^{S}_{a}$ and $\delta^{S}_{b}$ do not appear in Eqs. (3) since they contain a single $a_{n}$ term. Moreover, it is interesting to note that the terms proportional to $\delta^{S}_{a}$ correspond to the off-diagonal components of the wavevector-dependent effective permittivity in the approach proposed in Refs. Chebykin1 ; Chebykin2 , whereas the term proportional to $\delta^{S}_{b}$ is equivalent to a magnetic contribution and it can be strongly enhanced in the plasmonic-resonance condition ($|\epsilon_{||}|\gg 1$) Pendry_00 .

.3 Details on Numerical Simulations

In the results shown in Figs. 2 and 3, the rigorous Bloch-theory-based electromagnetic solution is implemented via a standard transfer-matrix approach Born , by solving an eigenvalue problem with phase-shift boundary conditions Joannopoulos .

The full-wave numerical results in Fig. 3 are obtained by means of the finite-element-based commercial software package COMSOL Multiphyics v. 3.4 Comsol . In this case, we consider a computational domain composed of a single unit-cell, terminated with periodic boundary conditions along the $x$-direction, and of finite size ($L=10$ $\mu$m, with vacuum regions of $0.2$ $\mu$m at both sides) and scattering boundary conditions at the entrance and the exit facets (orthogonal to the $z$-axis) of the domain. Moreover, we assume a TM plane-wave illumination (with $y$-directed magnetic field) normally impinging along the positive $z$-direction. We apply an adaptive discretization of the computational domain with a maximum mesh size of $0.1\lambda$ (which results into about $27000$ degrees of freedom), and utilize the UMFPACK direct solver, with default parameters. The field distribution in Fig. 2(c) is obtained by averaging the $y$-directed magnetic field over the unit cell. The spatial spectrum in Fig. 2(d) is subsequently obtained via a $512$-point fast-Fourier-transform.

Funding. U.S. Army International Technology Center Atlantic (Grant No. W911NF-14-1-0315).