Beyond local effective material properties for metamaterials

TL;DR

This paper investigates nonlocal effective material properties of metamaterials, emphasizing the importance of accounting for strong spatial dispersion to accurately describe their physical behavior and improve modeling accuracy.

Contribution

It introduces a general formulation for nonlocal constitutive relations in metamaterials, extending beyond local responses to better match actual dispersion properties.

Findings

Isofrequency surfaces match actual metamaterials better with nonlocal models

Strong spatial dispersion significantly affects metamaterial behavior

Nonlocal laws are necessary for accurate effective descriptions

Abstract

To discuss the properties of metamaterials on physical grounds and to consider them in applications, effective material parameters are usually introduced and assigned to a given metamaterial. In most cases, only weak spatial dispersion is considered. It allows to assign local material properties, i.e. a permittivity and a permeability. However, this turned out to be insufficient. To solve this problem, we study here the effective properties of metamaterials with constitutive relations beyond a local response and take strong spatial dispersion into account. The isofrequency surfaces of the dispersion relation are investigated and compared to those of an actual metamaterial. The significant improvement provides evidence for the necessity to use nonlocal material laws in the effective description of metamaterials. The general formulation we choose here renders our approach applicable to a wide class of metamaterials.

Tables

Table 1. Table 1: Parameter set of the SSD models for the best fit as shown in Fig.  5 . Depending on the model there is always one sign (the ± plus-or-minus \pm sign in Eqs. ( 15 ) and ( 26 )) that fits better to the fundamental Bloch mode. Bear in mind that parameters from different models have different expressions and dimensions.

Model
and sign:	4th(-)	SYM(+)
$p_0^{\text{TM}}$	$\left(-1.01-56.62i\right)\mu$m-2	$\left(-1.97-13.85i\right)\mu$m-2
$p_1^{\text{TM}}$	$\left(-2.35+0.28i\right)\cdot {10}^3\mu$m-2	$-7.21-1.12i$
$q_0^{\text{TM}}$	$0.804+0.807i$	$\left(-18.2-10.61i\right)\mu$m-2
$q_1^{\text{TM}}$	$-26.43-21.80i$	$\left(-27.7+33.95i\right)\mu$m-2

Equations

$$\mathbf{D}(\mathbf{r},\omega )$$

$$\mathbf{H}(\mathbf{r},\omega )$$

$$\mathbf{D}(\mathbf{r},\omega )={\int }_{{\mathbb{R}}^3}\hat{R}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )\mathbf{E}({\mathbf{r}}^{\prime },\omega )d{\mathbf{r}}^{\prime }$$

$$\mathbf{H}(\mathbf{r},\omega )={\mu }_0^{-1}\mathbf{B}(\mathbf{r},\omega ),$$

$$\mathbf{D}(\mathbf{r},\omega )={\int }_{{\mathbb{R}}^3}\hat{R}(\mathbf{r}-{\mathbf{r}}^{\prime },\omega )\mathbf{E}({\mathbf{r}}^{\prime },\omega )d{\mathbf{r}}^{\prime },$$

$$\tilde{\mathbf{D}}(\mathbf{k},\omega )=\hat{\tilde{R}}(\mathbf{k},\omega )\tilde{\mathbf{E}}(\mathbf{k},\omega ).$$

$${\tilde{D}}_{\alpha }(\mathbf{k},\omega )$$

$${\tilde{D}}_{\alpha }(\omega )=({\delta }_{\alpha \beta }+a_{\alpha \beta }){\tilde{E}}_{\beta },$$

$${ϵ}_{\alpha \beta }(\omega )={\delta }_{\alpha \beta }+a_{\alpha \beta }(\omega ).$$

$$D_{\alpha }(\mathbf{r},\omega )$$

$$b_{\alpha \beta \gamma }(\omega )\frac{\partial }{\partial x_{\gamma }}\stackrel{!}{=}[\hat{\xi }(\omega )\nabla \times {]}_{\alpha \beta }$$

$$c_{\alpha \beta \gamma \delta }(\omega )\frac{{\partial }^2}{\partial x_{\gamma }\partial x_{\delta }}\stackrel{!}{=}[\nabla \times \hat{\gamma }(\omega )\nabla \times {]}_{\alpha \beta }.$$

$$e_{\alpha \beta \gamma \delta \zeta \theta }(\omega )\frac{{\partial }^4}{\partial x_{\gamma }\partial x_{\delta }\partial x_{\zeta }\partial x_{\theta }}\stackrel{!}{=}[\nabla \times \nabla \times \hat{\eta }(\omega )\nabla \times \nabla \times {]}_{\alpha \beta },$$

$${\tilde{D}}_{\alpha }(\mathbf{k},\omega )=$$

$$+[\mathbf{k}\times \mathbf{k}\times \hat{\eta }\mathbf{k}\times \mathbf{k}\times {]}_{\alpha \beta }{\tilde{E}}_{\beta }.$$

$$\mathbf{k}\times \mathbf{k}\times \tilde{\mathbf{E}}+k_0^2\tilde{\mathbf{D}}[\tilde{\mathbf{E}}]=0.$$

$$k_x^2{\mu }_{xx}+k_z^2{\mu }_{zz}-{\left(k_x^2+k_z^2\right)}^2{\eta }_{yy}{\mu }_{xx}{\mu }_{zz}{k}_0^2={ϵ}_{yy}{\mu }_{xx}{\mu }_{zz}{k}_0^2,$$

$$k_z^2(k_x,k_0)$$

$$\pm \sqrt{{\left(p_0^{\text{TE }}\right)}^2-q_1^{\text{TE}}+2k_x^2\left(p_1^{\text{TE}}-p_0^{\text{TE}}\right)}$$

$$q_1^{\text{TE}}(k_0)$$

$$p_0^{\text{TE}}(k_0)$$

$$p_1^{\text{TE}}(k_0)$$

$$k_x^2{ϵ}_{xx}+k_z^2{ϵ}_{zz}-k_x^2{k}_z^2({ϵ}_{xx}{\eta }_{xx}+{ϵ}_{zz}{\eta }_{zz}){\mu }_{yy}{k}_0^2$$

$$-(k_x^4{ϵ}_{xx}{\eta }_{xx}+k_z^4{ϵ}_{zz}{\eta }_{zz}){\mu }_{yy}{k}_0^2={ϵ}_{xx}{ϵ}_{zz}{\mu }_{yy}{k}_0^2,$$

$$k_z^2(k_x,k_0)$$

$$\pm \sqrt{{\left[p_0^{\text{TM}}+\frac{q_0^{\text{TM}}-q_1^{\text{TM}}}{2}{k}_x^2\right]}^2-p_1^{\text{TM}}},$$

$$q_0^{\text{TM}}(k_0)$$

$$q_1^{\text{TM}}(k_0)$$

$$p_0^{\text{TM}}(k_0)$$

$$p_1^{\text{TM}}(k_0)$$

$$k_z^{{\text{WSD}}^2}(k_x,k_0)={\alpha }_1(k_0)+{\alpha }_2(k_0)k_x^2,$$

$${\alpha }_1^{\text{TE}}(k_0)$$

$${\alpha }_2^{\text{TE}}(k_0)$$

$${\alpha }_1^{\text{TM}}(k_0)$$

$${\alpha }_2^{\text{TM}}(k_0)$$

$$ϵ(x,y,z)=ϵ(\pm x,\pm y,\pm z).$$

Full text

I Introduction

Assigning effective material properties is a main problem in the macroscopic description of optical metamaterials Alù (2011); Smith (2010). The desire of replacing such materials by hypothetical, homogeneous ones rises from the simplification when discussing them within a physical framework. Once a metamaterial is homogenized, it might be described and discussed on the same level as a natural material by its effective material properties. Moreover, addressing numerically the wave propagation inside a homogeneous material is more efficient, i.e. less computational effort is required, than in performing a rigorous computation of the full structure. However, it is of utmost importance that both descriptions for the same metamaterial, i.e. the actual and the homogenized metamaterial, provide the same response, up to a certain accuracy, to the electromagnetic field. Otherwise the homogenization procedure is meaningless.

Various techniques have been established to assign effective material properties to optical metamaterials Andryieuski et al. (2012); Markel and Schotland (2012); Sihvola (2013); Silveirinha (2007); Liu and Alù (2013); Chipouline et al. (2012); Markel and Tsukerman (2013); Tsukerman and Markel (2016); Grahn et al. (2013a, b); Shevchenko et al. (2017). However, most previous techniques assume local constitutive relations, i.e. in a homogeneous material, the functional dependency of the auxiliary fields, ${\bm{D}}$ and ${\bm{H}}$, is a linear combination of the macroscopic fields ${\bm{E}}$ and ${\bm{B}}$, in which the coefficients are spatially independent Strunc (2007); Kong (2000)

[TABLE]

These equations are usually called local bi-anisotropic material laws, with effective material parameters $\hat{\epsilon}$, $\hat{\mu}$, $\hat{\xi}$ and $\hat{\zeta}$, where the hats refer to the quantities being tensors. A comprehensive review might be found in Ref. Semchenko et al., 2001. Although such constitutive relations are usually assumed, they exhibit several limitations Menzel et al. (2010); Zhukovsky et al. (2015); Andryieuski et al. (2015). In particular, the effective properties are only tensors that do not depend on the considered spatial frequency. However, they turn out to be inadequate when considering light propagation inside the metamaterial in an arbitrary direction and not just in the direction that was considered in the retrieval procedure Menzel et al. (2008, 2010). Moreover, whereas it can be safely expected that such constitutive relations are valid when considering metamaterials for which the operational wavelength is much longer than the size of the unit cell, they fail to be predictive for most metamaterials that are operated in a regime where the wavelength is not much smaller than size of the unit cell but only smaller. Such operational regime, unfortunately, is necessary to observe many relevant dispersive effects.

To overcome these limitations, we propose in this work two formulations for advanced constitutive relations in order to model metamaterials beyond a local response and to take strong spatial dispersion into account. The applicability of these models is investigated. The justification for such models is derived by studying here the isofrequency surfaces of the dispersion relation obtained by our models that are compared to those of an actual metamaterial when considering the fundamental mode, i.e. the mode with the smallest imaginary part in the propagation constant. We find clear evidence that it is necessary to consider strong spatial dispersion, i.e. nonlocal constitutive relations in the effective description of the metamaterials. We emphasise that in this contribution we study the bulk properties of the metamaterials in terms of the dispersion relation and/or isofrequency plots for the eigenmodes that are supported by the metamaterial. To make that approach viable requires also the formulation of additional boundary conditions to solve the problem how the modes couple at the interface from one medium to the other. These additional boundary conditions for the materials considered here are beyond the scope of this contribution and will be discussed elsewhere111A. Khrabustovskyi, K. Mnasri, M. Plum, C. Rockstuhl, C. Stohrer, “in preparation”. In the context of nonlocal constitutive relations, comparable approaches were already employed to effectively homogenize metamaterials beyond their local regime Ghasemi et al. (2015); Ciattoni and Rizza (2015); Fietz and Soukoulis (2012); Geng et al. (2015); Sun et al. (2015); Gorlach and Belov (2015); Silveirinha and Belov (2008); Tsukerman (2011). However, very often they required a specific geometry for the metamaterial that motivated the formulation of specific constitutive relations. In contrast, in this work we propose a scheme based on a phenomenological ansatz that can be applied to any periodic metamaterial made of subwavelength, resonant unit-cells to describe homogenized metamaterials.

This paper is structured as follows. In Sec. II we introduce our models from a very general ansatz. Taking a specific fishnet metamaterial as an example, we apply in Sec. III our models to describe its dispersion relation and to show the improvement, compared to a dispersion relation derived with a local materials law. Discussion and conclusion of this work are contained in Sec. IV.

II Constitutive Relations

Let us consider a periodic metamaterial in which the inclusions are intrinsically non-magnetic and reciprocal. We consider time-harmonic fields and only a linear response. It is therefore legitimate to assume that the electromagnetic response may be completely described by a non-local electric response, which can be written as the following integral form Landau et al. (1984)

[TABLE]

and

[TABLE]

where the kernel $\hat{R}$ represents the electric response tensor, that spatially links in a nonlocal fashion the electric field ${\bm{E}}$ to the displacement field ${\bm{D}}$. In a homogeneous medium, the response kernel does not depend on the exact spatial position but only on the difference between two considered points in space. This suggests that the kernel $\hat{R}$ in Eq. (2) reduces to a difference kernel, i.e.

[TABLE]

which is a convolution integral. It is more convenient to formulate this constitutive relation in spatial Fourier space, such that the convolution becomes a product

[TABLE]

This expression is too general for practical purposes. In order to reach useful constitutive relations, we assume that the unit cells are sub-wavelength and the kernel can be expanded into a Taylor expansion with respect to the spatial frequency ${\bm{k}}$. Up to the fourth order expansion, the kernel reads as the following

[TABLE]

where Einstein’s summation convention has been used. Note that all coefficients are in general complex functions of the frequency $\omega$. In the long wavelength limit, i.e. $|{\bm{k}}|\rightarrow 0$, spatial dispersion disappears and the constitutive relation (6) reduces to the one known from an anisotropic medium with only an electric response in its local description, i.e.

[TABLE]

with the linear electric polarization $\tilde{P}_{\alpha}=a_{\alpha\beta}\tilde{E}_{\beta}$. Therefore, the first line of Eq. (6) refers to the local permittivity of an anisotropic medium with

[TABLE]

In order to determine the nature of the higher order terms, we perform an inverse Fourier transform to real space. For plane waves, an inverse Fourier transform to the real space would lead to higher order derivatives of the electric field with respect to spatial coordinates at the same position where the displacement field is to be evaluated, whereupon the constitutive relation becomes

[TABLE]

where all coefficients were multiplied by a prefactor $(-\mathrm{i})^{n}$, with $n$ being the order of the spatial derivative. To be able to practically work with feasible constitutive relation, particular in the context of metamaterials, assumptions to these coefficients have to be made in order to reach well established effective medium theories that shall describe the actual metamaterial. For instance, a frequently made assumption is that

[TABLE]

and/or

[TABLE]

If these assumptions are met indeed and it is furthermore assumed that all higher order terms vanish, local constitutive relations identical to those in Eq. (1) can be derived by exploiting degrees of freedom that allow to suitably gauge Maxwell’s equations. The first assumption leads to a local optical activity (gyrotropy), whereas the second assumption leads to a local magnetic permeability Semchenko et al. (2001), respectively. We will refer to these assumptions as the weak spatial dispersion (WSD) approximation. “Spatial dispersion” because it results form non-local material laws and “weak” because there is a transformation that gives local constitutive relations, with a local magnetoelectric coupling and a local permeability, hence a local bi-anisotropic medium. However, it has been already shown that these assumptions are not enough to adequately describe the dispersion relation of an actual metamaterial beyond the paraxial regime Menzel et al. (2010). In a principal coordinates system, dispersion relations obtained from WSD are either of the hyperbolic kind, if the principal components of the permittivity or the permeability have opposite signs, or of the elliptic kind otherwise. However, most metamaterials with a nonlocal response usually show dispersion relations that widely differ from hyperbolas or ellipsoids. Exemplarily, we study in the following one example for a metamaterial where the dispersion relation of the fundamental mode indeed differs from that obtained within WSD. This suggests that the WSD approximation is not enough and we, therefore, need to go beyond it and have to retain higher order terms in the expansion. Here, we retain up to the fourth spatial-derivative of the electric field in Eq. (7).

In order to proceed, we assume WSD as a starting point and extend it into two directions. In our first model, we extend WSD by a, quite similar, but fourth order law in which the fourth order expansion coefficients take the following form

[TABLE]

where $\hat{\eta}=\text{diag}(\eta_{xx},\eta_{yy},\eta_{zz})$ is a diagonal matrix and represents a higher-order material parameter. This special form is chosen due to mathematical reasons. The coefficients that obey Eq. (10) do not change the variational class of Maxwell’s equations and hence, interface conditions may be found by means of variational methods. Even higher order terms again are assumed to vanish.

Our second model remains to be a second order law but it strictly orients on symmetry considerations of a unit-cell of a specific metamaterial. Specifically, it takes more coefficients into account than those that lead to local constitutive relations, i.e. assumptions (8) and (9). The model takes every coefficient into account as required by the symmetry of the considered metamaterial. The tensor elements of the kernel in Eq. (5) reflect the spatial symmetry of the actual structure. The degrees of freedom that the kernel might have can be reduced considerably depending on the symmetry of metamaterials. In this approach, we do not include fourth order terms in expansion (6), as it already gives significant improvements compared to WSD. The treatment of both approaches simultaneously, i.e. retaining both unit-cell’s symmetry conditions and a fourth order response, in principle, can be applied but the analysis becomes much more involved and is beyond the scope of this paper. In both models and for the sake of simplicity, we only consider metamaterials whose unit cells have a center of symmetry. By the presence of this (spatial inversion-)symmetry, all odd terms in expansion (7) vanish, hence no gyrotropic medium or optical activity is considered. However, it has to be stressed that this is not a general limitation. It is done here to concentrate only on one specific aspect.

II.1 Analysis by retaining fourth order response

Let us now investigate the dispersion relation and the consequences that arise from the following material law

[TABLE]

This is a direct extension of WSD by the next higher order response such that Maxwell’s equations remain of the varuatuinal class. Hence, boundary conditions may be found by means of variational methods.

In general, the dispersion relation can be derived by solving the wave-type equation

[TABLE]

This represents three coupled equations of second degree for the spatial frequency ${\bm{k}}$. Under the assumption that $k_{y}=0$ and considering the $z$-direction as the principle propagation direction, the dispersion relation for the transversal electric (TE) mode can be found. The dispersion relation expresses here the functional dependency of the frequency and the wave vector components. In TE-polarization the mode has a non-zero electric field component normal to the $k_{z}$-$k_{x}$-plane. The dispersion relation reads as

[TABLE]

where $\mu_{ii}=\frac{1}{1-k_{0}^{2}\gamma_{ii}}$. The solutions are

[TABLE]

with the frequency dependent coefficients

[TABLE]

where $q_{1}^{\text{TE}}$ has the dimension of m*-4* and both $p_{0}^{\text{TE}}$ and $p_{1}^{\text{TE}}$ the dimension of m*-2*. For the transversal magnetic (TM) mode, i.e. the mode that has an electric field in the $k_{x}$-$k_{z}$-plane, we obtain

[TABLE]

with the solutions

[TABLE]

where

[TABLE]

Here, both $q_{0}^{\text{TM}}$ and $q_{1}^{\text{TM}}$ are dimensionless, while $p_{0}^{\text{TM}}$ and $p_{1}^{\text{TM}}$ have the dimensions of m*-2* and m*-4*, respectively. Due to the higher number of degrees of freedom, the functional dependency of $k_{z}(k_{x},k_{0})$ is more advanced, in comparison to the previously proposed relation Menzel et al. (2010),

[TABLE]

which represents the dispersion relation of a material with local constitutive relations (WSD). Here, the coefficients $\alpha_{1}$ and $\alpha_{2}$ depend on the polarization where for the TE mode they read

[TABLE]

and

[TABLE]

for the TM mode. In general, isofrequency contours of media described by Eq. (17) are limited to two cases: they are either of hyperbolic or elliptic kind. This limitation is lifted by introducing more complicated dispersion relations, e.g., our fourth order response. To illustrate the possible features an isofrequency contour may have when considering such fourth order constitutive relations, Fig. 1 shows some plots of $k_{z}(k_{x})$ for a fixed frequency $k_{0}$ for some generically chosen parameters. The obtained isofrequency contours give rise to more advanced curves. They allow for a homogenous description of an metamaterial with disperive features inaccessible by a local material.

II.2 Analysis with respect to structure’s symmetry

In this model, instead of taking higher orders in the expansion of the kernel into account, we take a deeper look in the geometry of an actual (real) structure. As an example we consider the fishnet metamaterial (see Fig. 3) as a subject of homogenization. The unit cell of the fishnet metamaterial is symmetric under transformation with respect to three orthogonal mirror planes, i.e. the permittivity distribution obeys

[TABLE]

This symmetry class is also known as orthorhombic symmetry Agranovich and Ginzburg (2014), noted as D${}_{2\text{h}}$. According to this consideration, we can write the expansion of the displacement field in the form

[TABLE]

where Einstein’s summation rule applies only to $\beta$. Here, we assume that the local permittivity and $\hat{\gamma}$ are diagonal. More importantly, some of the second order terms may be written in a similar way as in the WSD. The only difference to the WSD is the higher order susceptibility contribution $c_{\alpha\alpha\alpha\alpha}$ that couples $k_{\alpha}^{2}$ with $E_{\alpha}$. This term, however, cannot be set to zero a priori for a metamaterial with the considered symmetry of the fishnet. It has to be taken into account and appears on the same footing as the other second order terms. These other terms in $k^{2}$ are assumed to obey the condition that yields a local magnetic response, i.e. Eq. (9). The dispersion relation follows from solving the wave equation as described in Eq. (12). We assume equally that $k_{y}=0$ and we study a plane with a wavevector in the $k_{x}$-$k_{z}$-plane. The dispersion relation for the TE mode is then

[TABLE]

with the solutions

[TABLE]

where

[TABLE]

The fact that the TE mode does not experience any strong spatial dispersion relies on the nature of the expansion in Eq. (LABEL:eq:constitD2H) with the nonlocal response that couples only in the direction of the electric field, hence no cross coupling between the displacement field and the electric field as can be seen in Eq. (LABEL:eq:constitD2H). In contrast, a more complicated dispersion relation is found for the TM modes with

[TABLE]

where $\mu_{ii}=\frac{1}{1-k_{0}^{2}\gamma_{ii}}$ and $c_{J}=c_{jjjj}$. This equation is biquadratic in $k_{x}$ and $k_{z}$ and describes more complicated isofrequencies than quadratic equations, e.g. given by WSD. Of course, the limit $(c_{X},c_{Z})\rightarrow(0,0)$ restores the dispersion relation given by WSD. The solutions of Eq. (25) are

[TABLE]

with the related coefficients that read

[TABLE]

It has to be noted that all coefficients have the dimension of m*-2*, except $p_{1}$ being dimensionless. Here, we have four independent coefficients which increase the number of degrees of freedom to four. Previously, e.g. in Eq. (17), the dispersion relation contained only two independent coefficients, $\alpha_{1}$ and $\alpha_{2}$, hence only two degrees of freedom. Moreover, and this holds for both Strong Spatial Dispersion (SSD) models, each of the Eqs. (26, 13, 15) yield four possible solutions for $k_{z}(k_{x},k_{0})$ from which two with positive and two with negative imaginary parts. We consider here only solutions with a positive imaginary part as they describe exponentially damped solutions in our principle propagation direction. In contrast the two solutions with a negative imaginary part would correspond to exponentially growing solutions. They are unphysical for a passive medium. However, this still suggests that in the actual homogeneous medium characterized by material laws beyond WSD, more than a single mode is excited at the interface where the continuity of the tangential wave vector components dictates which modes are excited. To simplify the analysis, we concentrate on investigating the fundamental mode, i.e. the mode with the smallest positive imaginary part. In general, the imaginary parts of the two solutions with $\Im(k_{z})>0$ may cross and a mode transition has to be taken into account. Here as well, Fig. 2 shows some isofrequency contours from selected parameter sets. The complexity of their shapes suggests the ability to capture the effects of strong spatial dispersion.

In the next section, we show the importance of retaining these nonlocal effects in the effective description of the metamaterial by directly showing the improvements that follow from taking SSD into account, as introduced in Eqs. (11) and (LABEL:eq:constitD2H).

III Application to a fishnet metamaterial

In this section, a numerical experiment is done to get access to the dispersion relation of an actual structure as a reference. We consider the fishnet metamaterial shown in Fig. 3 as an example, which is known to exhibit a negative refractive index at optical frequencies. The geometrical parameters are taken from literature Dolling et al. (2006). The unit cell’s dimensions are $\Lambda_{x}=\Lambda_{y}=600$ nm, $\Lambda_{z}=200$ nm. The rectangular holes are made of perpendicularly aligned nanowires with thicknesses of $w_{x}=100\,$nm and $w_{y}=316\,$nm. It is made of two $45\,$nm Silver (Ag) layers where their permittivity obey the Drude model for metals, with the plasma frequency being $\omega_{p}=13700\text{ THz}$ and the relaxation rate $\Gamma=85\text{ THz}$. These layers are separated by a $30\,$nm Magnesium fluoride (MgF2) spacer whose permittivity is assumed to be nondispersive ($\epsilon_{\text{MgF}_{2}}=1.9044$) in the frequency range of interest. Furthermore, the unit cell is symmetric with respect to spatial inversion, i.e. $\epsilon({\bm{r}})=\epsilon(-{\bm{r}})$.

In order to calculate the Bloch mode dispersion relation, i.e. $k_{z}=k_{z}(k_{x},k_{y},k_{0})$ where $k_{0}$ represents the wave number in free space, a plane-wave expansion ansatz is numerically performed. In general, $k_{z}$ can be complex. Its real part, $\Re(k_{z})$, refers to the oscillatory part while its imaginary part, $\Im(k_{z})$, denotes the energy loss in the principal propagation direction. We restrict our considerations to the solutions with $\Im(k_{z})>0$. The dominating Bloch mode, i.e. the fundamental Bloch mode that prevails after a finite propagation length is the one with the smallest positive $\Im(k_{z})$ as all higher modes experience much stronger damping. Fig. 4 shows both dispersion relation and isofrequency contours of the fundamental Bloch mode for different transverse vectors $k_{x}\geq 0$ at a fixed frequency and for different frequencies at a fixed transverse wave vector, respectively. Here we restrict ourselves to the $k_{x}$-$k_{z}$-plane, i.e. $k_{y}=0$. The mode is TM polarized. For wavenumbers around $4.3\,\mu$m*-1*, we observe negative $\Re(k_{z})$, which implies a negative index, i.e. momentum and energy flux propagate in opposite directions. We are basically interested in homogenizing the material around this resonance wavenumber. Figure 4 (right) shows the isofrequency contour for $k_{0}=4.3\,\mu$m*-1* of the real structure. This numerically obtained isofrequency contour of the actual metamaterial has to be reproduced by the dispersion relation of the effective medium from WSD as well as from both SSD with some fixed set of parameters. The comparison is based on a least absolute deviations fit by optimizing the parameters $(q_{0},p_{0},q_{1},p_{1})$ of fundamental TM modes from the SSD models and $(\alpha_{1},\alpha_{2})$ from the WSD. As a quantity of measure, we define the merit function as the following

[TABLE]

with a suitably defined weight function $w(k_{x})$. Here $\text{i}\in\{\text{WSD, 4th, SYM}\}$ denotes the model taken into account. 4th and SYM refer to the fourth-order and to the symmetry model, respectively. As all the expressions were derived from a Taylor expansion for small ${\bm{k}}$ (see Eq. (6)), it is legitimate to introduce a weight $w(k_{x})$ such that the fitting procedure is more focused for small $k_{x}$. Here we chose an exponentially decreasing dependency, i.e.

[TABLE]

where $\alpha$ was chosen to be $\alpha=2.5\Lambda_{x}$, with $\Lambda_{x}=0.6\,\mu$m being the lateral period of the fishnet structure. The results for a selected frequency – here we chose the worst case scenario, i.e. at the resonance – when considering the optimized parameters are depicted in Fig. 5. It shows the isofrequency contour in both real and imaginary part of the dispersion relation numerically calculated for the actual fishnet metamaterial and the dispersion relation obtained for the best fit of parameters at the resonance frequency for the three different models considered. The parameter set of the SSD models for the best fit and the right signs are summarized in Table 1. For the parameters of WSD, i.e relation (17), we obtain $\alpha_{1}^{\text{TM}}=(-5.85-28.12\mathrm{i})\mu$m*-2* and $\alpha_{2}^{\text{TM}}=-2.49+3.79\mathrm{i}$. Clearly, WSD is only in a good agreement with the isofrequency contour of the real structure in the paraxial regime, i.e. $k_{x}\ll k_{0}$. Beyond the paraxial regime, we recognize from the shape of the black curves of $k_{z}(k_{x})$, that WSD, which by nature is either an ellipse or a hyperbola, is not enough to describe the dispersion relation of such complicated structure. This limitation can be lifted by considering nonlocal constitutive relations as proposed above. The actual dispersion relation can be much better described when homogenizing the metamaterial with constitutive relations beyond the local model. To quantify the actual improvement, we study the merit function as a function of the frequency for the three different constitutive relations obtained in here. Each value for the merit function has been obtained from an individual fit at a specific frequency.

Figure 6 shows the improvements in effectively describing the metamaterial with nonlocal constitutive relations for all of the simulated frequencies. Both SSD models are more accurate than WSD. The integrated error that expresses how good a specific constitutive relation in a homogenized medium can explain the actual dispersion relation of the given metamaterial, is in average two orders of magnitude better for the nonlocal material laws. In resonance, the deviation is strongest irrespective of the considered model. However, this is somehow expected that the effective description tends to be inaccurate in the resonance regime. Nevertheless, the findings immediately implies retaining nonlocal constitutive relations is required for a more realistic homogenization of metamaterials.

IV Discussion and Conclusions

The aim of this paper was to introduce on the one hand a viable route to describe homogenized metamaterial beyond the assumption of a local response and to show the importance of retaining nonlocal material laws in order to adequately describe the dispersion relation of a homogenized metamaterial. Of course we are not the first to work in this direction, but the approach we have choosen here is generally applicable and does not hinge on the assumption of a specific geometry of the metamaterial. At the example of a fishnet metamaterial we have shown that WSD is not enough to properly capture the dispersion relation as their functional dependency and their isofrequency contours – either hyperbolic or elliptic – are too simplistic to give accurate predictions beyond the paraxial regime. Significant improvement only comes by introducing nonlocal material laws in their effective description. These come with more degrees of freedom, hence more complicated isofrequency come into play and give a better description of metamaterials. In our work we have studied the light propagation in the bulk of a metamaterial and we obtained all the coefficients that enter the dispersion relation by means of optimization in comparision to the numerically calculated dispersion relation of an actual metamaterial. However, for a complete inversion, i.e. to obtain the individual material parameters, requires finding the interface conditions between one ordinary medium and one with nonlocal material laws. This is subject of ongoing research.

Acknowledgements

We gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) through CRC 1173. K.M. also acknowledges support from the Karlsruhe School of Optics and Photonics (KSOP). We would like to thank Radius Suryadharma and Andreas Vetter for proofreading the manuscript.