Generalized method for retrieving effective parameters of anisotropic metamaterials

A. Castanié; J.-F. Mercier; S. Félix; A. Maurel

doi:10.1364/OE.22.029937

1. Introduction

Metamaterials have generated an intense research interest in recent years for their ability to produce unusual properties of wave propagations, in electromagnetism, acoustics, elasticity and for water waves, see, e.g., [1–5]. These structures are most commonly made of a periodic arrangement of unit cells, that may be very complex, and there is a need for proper characterization of the resulting behavior at large scale in order to better understand and further exploit them.

The notion of effective parameters for metamaterials has been introduced in analogy with the notion of homogeneous parameters defined for natural materials, being inhomogeneous at small scales. This is a simplifying approximation expected to be valid when the scales of the material inhomogenities are much smaller than the scale associated to the spatial variation in the acoustic, or electromagnetic, field. In the low frequency regime, homogenization theory provides the effective parameters unambiguously (see, e.g, [6]) and recent studies have shown the possibility to extend such approach to higher frequencies, typically when the unit cell presents resonances [7]. Alternatively, and we may say more pragmatically, a popular method consists in determining the effective parameters using an inversion technique [8–11] assuming that the reflection and transmission coefficients by a metamaterial slab illuminated by a plane wave are known. The idea beyond such methods is to get material properties with broadband validity and to avoid mathematical approaches hardly suitable in practice. For resonant structures, this means that the dependence of the effective parameters on the frequency are looked for, and this may lead to negative effective parameters [11–15]. For non resonant structures, the parameters are expected to remain relatively unchanged but the idea is to get more accurate parameter values in the frequency range of interest [16–18] (more accurate than the values obtained from homogenization). As pointed out in [19], the main virtue of this technique is that it does not require any assumption except that the concept of effective medium is valid, which is not incidental [20–24].

The retrieval method leads to the effective refractive index n and the impedance ratio ξ (between, say, air and the effective medium) using the reflection and transmission coefficients through a thin slab containing few unit cells of the metamaterial. In its initial form, the method was adapted to isotropic effective medium and the method used the scattering properties of a slab illuminated by plane waves incident normally to the slab interfaces. Firstly proposed by Smith in 2002 [12] in the context of electromagnetic waves, the method has been improved to the case of an asymmetric slab [9], to the case of an isotropic effective medium with obliquely incidences [25] and to the case of an anisotropic effective medium with normal incidence [26] and with oblique incidence [27]; however, in [26, 27], it was assumed that the interfaces were parallel to one principal axis of anisotropy. In addition to these improvements, the ambiguities in the inversion technique have been clearly addressed [19, 28].

Note that in 2007, Fokin and coauthors propose a slightly different retrieval method [11] in the context of acoustic waves. Although the situation in both contexts of waves is strictly the same in two dimensions (except there is no polarization option in acoustics), the two literatures, in acoustic and electromagnetism, remain largely separated.

In this paper, the main remaining restriction of the retrieval methods, namely the assumption that the principal directions of anisotropy coincide with the directions of the unit cell, is addressed (angle α in Fig. 1). This is done together considering oblique incidence of the plane wave on the air-metamaterial interface (angle θ). This may be of practical importance. Indeed, retrieval methods have been applied to characterize anisotropic metamaterial with α = 0, and have revealed interesting phenomena: the design of anisotropic metamaterials with negative index operating at variety of incidence angles [29], the highly directive emission of a source embedded in an anisotropic metamaterial [30, 31], the design of a slab lens with low loss using an anisotropic metamaterial with negative properties involving pure dielectric rod-type structure, the rainbow-like emission of a source placed in a uniaxial metamaterial slab [32], the inversion of the Brewster angle in anisotropic metamaterials with negative properties. We think that the ability to inspect effective properties when the resulting anisotropy is not aligned with the unit cell (thus, α ≠ 0) may open doors to interesting applications.

Fig. 1 (a) Typical configuration of the real medium, the reflection R and transmission T through a slab containing few cells (otherwise ℓ_x-periodic in the x-direction) are calculated for an obliquely incident plane wave with angle θ and wavenumber k. (b) Effective anisotropic medium leading to the same scattering properties (R and T).

Download Full Size | PDF

The paper is organized as follows. In Sec. 2, the 2D wave equations for electromagnetism and acoustic are briefly recall, to clarify the equivalences. In Sec. 3, the retrieval method for anisotropic effective media and oblique incidence is presented. This is done following [11] and by considering first the case α = 0 (Sec. 3.1). Although a similar calculation has been already done in [27], this “warm up” allows to better understand the case α ≠ 0 (presented in Sec. 3.2). Indeed, contrary to what is said in [28], the case of an isotropic slab (even at oblique incidence) has strong differences with the case of an anisotropic slab with arbitrary principal axes. Sec. 4 presents results in the particular case of layered media, for which available results coming from the homogenization theory have been tested [33]. In this particular configuration, our retrieval technique is validated, and the variations of the retrieved effective parameters with the frequency and the anisotropy angle are inspected. We also report in Sec. 5 results of the retrieval method for array of tilted elliptical particles (the acoustic case is considered here).

2. 2D Wave equation in anisotropic media for electromagnetic and acoustic waves

In this section, the 2D wave equations for wave propagating in an anisotropic medium is briefly derived, for electromagnetic polarized waves and for acoustic waves. A time harmonic dependence e^−iωt is assumed in the following.

First, we define in electromagnetism the permittivity tensor ε^diag and permeability tensor μ^diag, and in acoustics, the mass density tensor ρ^diag (the bulk modulus B is assumed to be scalar) expressed in the principal directions (X, Y, Z)

\begin{matrix} ε^{diag} = (\begin{matrix} ε_{X} & 0 & 0 \\ 0 & ε_{Y} & 0 \\ 0 & 0 & ε_{Z} \end{matrix}), & μ^{diag} = (\begin{matrix} μ_{X} & 0 & 0 \\ 0 & μ_{Y} & 0 \\ 0 & 0 & μ_{Z} \end{matrix}), & in electromagnetism, \\ ρ^{diag} = (\begin{matrix} ρ_{X} & 0 & 0 \\ 0 & ρ_{Y} & 0 \\ 0 & 0 & ρ_{Z} \end{matrix}), & in acoustics . \end{matrix}

The equations can be solved in a (x, y, z = Z) space deduced from (X, Y, Z) by a rotation R of angle α in the (x, y) plane,

R = (\begin{matrix} cos α & sin α & 0 \\ - sin α & cos α & 0 \\ 0 & 0 & 1 \end{matrix})

and we get

\begin{array}{l} in electromagnetism {\begin{array}{l} \nabla \times E = i k ({}^{t}R μ^{diag} R) H, \\ \nabla \times H = - i k ({}^{t}R μ^{diag} R) E, \end{array} \\ in acoustics {\begin{array}{l} \nabla p = i k ({}^{t}R ρ^{diag} R) u, \\ B \nabla \cdot u = i k p, \end{array} \end{array}

which correspond to the Maxwell equations for the electric field E and the magnetic field H, and to the Euler equations for the pressure p and velocity u, respectively. In the 2D case, the medium does not couple the electric E and magnetic H fields and for polarized waves, either Transverse Magnetic (TM) for which H = (0, 0, H) or Transverse Electric (TE) for which E = (0, 0, E). A wave equation for respectively E or H similar to Eq. (5) is obtained, namely,

\begin{array}{l} \nabla \cdot [(\begin{array}{l} 1 / ε_{x} & 1 / ε_{x y} \\ 1 / ε_{x y} & 1 / ε_{y} \end{array}) \nabla H] + μ_{Z} k^{2} H = 0, & TM waves, \\ \nabla \cdot [(\begin{array}{l} 1 / μ_{x} & 1 / μ_{x y} \\ 1 / μ_{x y} & 1 / μ_{y} \end{array}) \nabla E] + ε_{Z} k^{2} E = 0, & TE waves . \end{array}

In the context of acoustics, 2D case corresponds to a configuration where the pressure field p(x, y) and velocity u(x, y) do not depend on z and we get

\nabla \cdot [(\begin{array}{l} 1 / ρ_{x} & 1 / ρ_{x y} \\ 1 / ρ_{x} & 1 / ρ_{y} \end{array}) \nabla p] + \frac{k^{2}}{B} p = 0,

with the definitions and correspondances indicated in the Table 1.

Table 1. Coefficients in the 2D wave equation for anisotropic media in electromagnetism and acoustics.

View Table

3. Retrieval method

The basic idea of the retrieval method lays in the assumption that the scattering properties of the medium, Fig. 1(a), can be described by the scattering properties of an equivalent effective homogeneous anisotropic medium, Fig. 1(b).

We start with the warm-up case α = 0, already considered in [27]. First, (R, T) are determined in the configuration of Fig. 1(b), afterwards the retrieval method is presented. Then, the case α ≠ 0 is presented.

3.1. Case of anisotropy direction α = 0

To begin with, the retrieval method is presented assuming α = 0. Outside the slab, say, in air, the wave field H satisfies (Δ + k²)H = 0. In the slab, that contains several unit cells in [0, ℓ_y], the wavefield satisfies

\nabla \cdot [(\begin{array}{l} 1 / ε_{Y} & 0 \\ 0 & 1 / ε_{X} \end{array}) \nabla H] + μ_{Z} k^{2} H = 0,

and (ε_X, ε_Y, μ_Z) are the permittivities and permeability relative to the permittivity ε^ext and permeability μ^ext outside the slab (and k² = ω²ε^ext μ^ext). The wave field inside and outside the slab can be written

{\begin{array}{l} H (y \leq 0) = e^{i k sin θ x} [e^{i k cos θ y} + R e^{- i k cos θ y}], \\ H (0 < y \leq ℓ_{y}) = e^{i k sin θ x} [a e^{i n k y} + b e^{- i n k y}], \\ H (ℓ_{y} \leq y) = e^{i k sin θ x} T e^{i k cos θ (y - ℓ_{y})} \end{array}

where nk = k_Y denotes the horizontal effective wavenumber (n is the index of the effective medium). Applying the continuity, at each interfaces y = 0 and y = ℓ_y, of the field H and ∂_yH/ε_X (inside the slab) and ∂_yH (outside the slab), we get the reflection R and transmission T coefficients

R = \frac{(1 - ξ^{2}) (e^{i n k ℓ_{y}} - e^{- i n k ℓ_{y}})}{{(1 + ξ)}^{2} e^{- i n k ℓ_{y}} - {(1 - ξ)}^{2} e^{i n k ℓ_{y}}}, T = \frac{4 ξ}{{(1 + ξ)}^{2} e^{- i n k ℓ_{y}} - {(1 - ξ)}^{2} e^{i n k ℓ_{y}}}

where R = R(θ), T = T (θ) and, for TM waves,

n \equiv \frac{k_{Y}}{k} = \sqrt{μ_{Z} ε_{X} - \frac{ε_{X}}{ε_{Y}} {sin}^{2} θ}, ξ = \frac{ε_{X} cos θ}{n} .

In [34], the case of TE waves is considered. The authors obtain (with our notations)

n = \sqrt{ε_{Z} μ_{X} - \frac{μ_{X}}{μ_{Y}} {sin}^{2} θ}, and ξ = \frac{μ_{X} cos θ}{n},

and ξ is used to calculate the reflection coefficient for a single interface R_hs = (ξ − 1)/(ξ + 1), afterwards R is deduced (see Eq. (3) of [34]).

At this stage, let us just make a parenthesis and stress the analogy with the acoustics case:

n = \sqrt{\frac{ρ_{Y}}{B} - \frac{ρ_{Y}}{ρ_{X}} {sin}^{2} θ}, and ξ \equiv \frac{ρ_{Y} cos θ}{n} .

The above expressions are in agreement with the equivalents given in the Table 1.

Next, we come back to the retrieval method. Following [11], it is easy to show that n and ξ can be calculated from (R, T) owing to the relations

{\begin{array}{l} n = \frac{i}{k ℓ_{y}} [log \frac{(1 - R^{2} + T^{2} + \tilde{ξ})}{2 T} + 2 i m π], with \tilde{ξ} \equiv \pm \sqrt{{(1 - R^{2} + T^{2})}^{2} - 4 T^{2}}, \\ ξ = \frac{\tilde{ξ}}{{(1 - R)}^{2} - T^{2}}, \end{array}

with m an integer. The determination of the sign in ξ̃ is known as “sign ambiguity”, and the problem of the 2π-multivaluation of n (choice of m) as “branch ambiguities”. These ambiguities have been addressed in [19,28] and they are disregarded here. For weak scattering, these ambiguities are easy to solve. Indeed, we expect in that case R ∼ 0 and T ∼ e^ikcosθℓ_y, associated to an almost perfect impedance matching ξ ∼ 1 and to almost equal horizontal wavenumbers: nk in the slab and k cosθ outside, namely n ∼ cosθ. Then, ξ̃ = [(1 − R)² − T²]ξ ∼ (1 − e^{2ikcosθℓ_y}) (which answers the sign ambiguity) leading to n ∼ cosθ in Eq. (12) if the principal determination is considered, (branch m = 0).

The steps of the inversion towards the effective parameters are as follows:

n(θ) and ξ (θ) are computed from (R, T) according to Eq. (12).
ε_X = ξ (θ)n(θ)/cosθ, from Eq. (9), is deduced.
n and ε_X being known, the function n²/ε_X = C₁ − C₂ sin² θ can be fitted as a function of θ to get μ_Z = C₁ and ε_Y = 1/C₂, Eq. (9).

In the following, we adapt this inversion to the case where α ≠ 0.

3.2. Case of tilted principal axes, α ≠ 0

If α ≠ 0, the wave field H satisfies the wave equation (Eq. (3))

\nabla \cdot [(\begin{array}{l} 1 / ε_{y} & 1 / ε_{x y} \\ 1 / ε_{x y} & 1 / ε_{x} \end{array}) \nabla H] + μ_{z} k^{2} H = 0,

with

{\begin{array}{l} \frac{1}{ε_{x}} = \frac{{cos}^{2} α}{ε_{Y}} + \frac{{sin}^{2} α}{ε_{X}}, \\ \frac{1}{ε_{y}} = \frac{{sin}^{2} α}{ε_{Y}} + \frac{{cos}^{2} α}{ε_{X}}, \\ \frac{1}{ε_{x y}} = - sin α cos α (\frac{1}{ε_{X}} - \frac{1}{ε_{Y}}) . \end{array}

The main difference with the previous calculation of R and T is that the horizontal wavenumber of the right- and left- going waves in the slab are not just of opposite sign. Denoting them k^±, the field inside the slab can be written

H (0 < y \leq ℓ_{y}) = e^{i k sin θ x} [a e^{i k^{+} y} + b e^{i k^{-} y}] .

Inserting this form in the wave equation (13), we get the dispersion relation

{k^{\pm}}^{2} / ε_{y} + 2 k^{\pm} k sin θ / ε_{x y} + k^{2} ({sin}^{2} θ / ε_{x} - μ_{Z}) = 0,

leading to

k^{\pm} = - \frac{ε_{y}}{ε_{x y}} k sin θ \pm n k, with n \equiv \sqrt{μ_{Z} ε_{y} - \frac{ε_{y}^{2}}{ε_{X} ε_{Y}} {sin}^{2} θ},

where we have used

1 / (ε_{x} ε_{y}) - 1 / ε_{x y}^{2} = 1 / (ε_{X} ε_{Y})

. At the interfaces y = 0 and y = ℓ_y, we apply the continuity of the field H and the continuity of ∂_yH outside the slab and ∂_yH/ε_y + ∂_xH/ε_xy inside the slab. We get

R = \frac{(1 - ξ^{2}) (e^{i n k ℓ_{y}} - e^{- i n k ℓ_{y}})}{{(1 + ξ)}^{2} e^{- i n k ℓ_{y}} - {(1 - ξ)}^{2} e^{i n k ℓ_{y}}}, T = \frac{4 ξ e^{- i k ℓ_{y} sin θ ε_{y} / ε_{x y}}}{{(1 + ξ)}^{2} e^{- i n k ℓ_{y}} - {(1 - ξ)}^{2} e^{i n k ℓ_{y}}}

and in this case, n and ξ are defined by

n \equiv \sqrt{μ_{Z} ε_{y} - \frac{ε_{y}^{2}}{ε_{X} ε_{Y}} {sin}^{2} θ}, ξ \equiv \frac{ε_{y} cos θ}{n} .

Obviously, α = 0 leads to 1/ε_xy = 0, ε_x = ε_Y and ε_y = ε_X, which is consistent with the previous calculation.

To extend the retrieval method to this case, the main problem is that ε_xy appears only in the expression of T through a phase term that is unknown a priori. Thus, it will not be possible to build combinations of R and T able to isolate n or ξ only, as it is done for α = 0, Eq. (9). This phase term has to be eliminated first, before any combination of R and T. This is a preliminary step in the inversion, which consists in building a phase compensated transmission coefficient T_c, and this requires T (θ) and T (−θ) to be known. If so, it is sufficient to remark that θ → −θ affects the extra phase term but does not affect ξ and n, and this leads to

e^{- 2 i k ℓ_{y} sin θ ε_{y} / ε_{x y}} = \frac{T (θ)}{T (- θ)},

from which we deduce

T_{c} (θ) \equiv \sqrt{T (- θ) T (θ)} = \frac{4 ξ}{{(1 + ξ)}^{2} e^{- i n k ℓ_{y}} - {(1 - ξ)}^{2} e^{i n k ℓ_{y}}},

Now, the inversion procedure can be apply straightforwardly to (R, T_c) following the same steps as described in Sec. 3.1:

n(θ) and ξ (θ) are computed from (R, T_c) according to Eq. (12)
ε_y = ξ (θ)n(θ)/cosθ, from Eq. (19), is deduced.
n and ε_y being known, the function n²/ε_y = C₁ − C₂ sin² θ can be fitted as a function of θ to get μ_Z = C₁ and ε_X ε_Y = ε_y/C₂.

At that stage, ε_y and the product ε_X ε_Y have been determined and two more steps are needed to get (ε_xy,ε_x):

From Eq. (20), ε_xy can be calculated, and this may be done by fitting Eq. (20) or by using a direct inversion (ε_xy = −2ikℓ_y sinθε_y/log[T (θ)/T(−θ)]). Note that a fitting procedure is possible here since ε_xy is expected to be independent of θ (at least in a certain range of θ), which is not the case for n in Eq. (12).
Finally, ε_x is determined using $\frac{1}{ε_{x}} = \frac{ε_{y}}{ε_{X} ε_{Y}} + \frac{ε_{y}}{ε_{x y}^{2}} .$ The full tensor can be described equivalently by (ε_X, ε_Y, α) owing to ${\begin{cases} \frac{1}{ε_{X}} = \frac{1}{2} [\frac{1}{ε_{x}} + \frac{1}{ε_{y}} + \sqrt{{(\frac{1}{ε_{x}} - \frac{1}{ε_{y}})}^{2} + \frac{4}{ε_{x y}^{2}}}], \\ \frac{1}{ε_{Y}} = \frac{1}{2} [\frac{1}{ε_{x}} + \frac{1}{ε_{y}} - \sqrt{{(\frac{1}{ε_{x}} - \frac{1}{ε_{y}})}^{2} + \frac{4}{ε_{x y}^{2}}}], \end{cases}$ where we assume ε_X ≤ ε_Y and ${\begin{array}{l} α = {sin}^{- 1} [\sqrt{\frac{1 / ε_{Y} - 1 / ε_{x}}{1 / ε_{Y} - 1 / ε_{X}}}], & if ε_{x y} \leq 0 \\ α = 180 ° - {sin}^{- 1} [\sqrt{\frac{1 / ε_{Y} - 1 / ε_{x}}{1 / ε_{Y} - 1 / ε_{X}}}], & if ε_{x y} > 0, \end{array}$ where the determination of α is done in [0, 180°]. The choice to assign the largest permittivity to the axis X in Eq. (23) is not restrictive since it only reflects the equivalence between the situations, for say ε₁ < ε₂, (ε_X = ε₁, ε_Y = ε₂, α) and (ε_Y = ε₁, ε_X = ε₂, 90° + α).

In the following sections, we illustrate our retrieval method in two cases: (1) the case of layered medium and we inspect the influence of the anisotropy direction α and of the wavenumber k on the retrieved parameters (electromagnetic wave are considered) and (2) the case of arrays of tilted elliptical particles and we inspect the influence of the array size, beyond the dominant effects of the filling fraction and of the orientation of the particles (acoustic waves are considered).

4. Results for layered dielectric media

The results presented in this section corresponds to a structure made of layers (slab, see Fig. 2). The layers are composed of alternating air and a dielectric material chosen with parameters relative to those of air ε₁ = 10 and μ₁ = 0.2, with periodicity d and filling fraction φ = 0.5. In the low frequency regime, layered structures are in general well described by homogenization theory [6, 33]; for comparison with the retrieved effective parameters, we give the values of the homogenized parameters

{\begin{array}{l} ε_{X}^{homog} = {(1 - φ + φ / ε_{1})}^{- 1} = 1.818, & ε_{Y}^{homog} = 1 - φ + φ ε_{1} = 5.500, \\ μ^{homog} = 1 - φ + φ μ_{1} = 0.6, \end{array}

with the notation of Eq. (6) (for α = 0, ε_x = ε_Y and ε_y = ε_X). To apply the retrieval method, one consider a slab composed of two unit cells (ℓ_x = d/cosα, ℓ_y = 2d/sinα when α ≠ 0 -see Fig. 2 - and ℓ_x = d, ℓ_y = 2d when α = 0). Along the x-direction, periodic, or Floquet, boundary conditions are imposed. Numerical calculations are performed using a multimodal method, similar to the Rigorous Coupled Wave Approach [35, 36]. From the numerics, the (R, T) coefficients are calculated, and allow then to deduce the retrieved parameters.

Fig. 2 Configuration of the study. The slab contains slanted layers, made of alternating air and a material with parameters relative to those of the air ε₁ = 10 and μ₁ = 0.2; the wave is incident with angle θ and we measure R and T.

Download Full Size | PDF

4.1. Validation of the retrieval method

To begin with, and to anticipate on the validation of our retrieval method, Figs. 3 and 4 show the comparisons between two wavefields, in the cases α = 0 and α = 45°, respectively: the wavefield computed considering the real structure (direct numerics, Figs. 3(a) and 4(a)), and the wavefield computed in an effective anisotropic slab of length L with parameters deduced from the retrieval method (Figs. 3(b) and 4(b)). For α = 0 and kd = 0.5, we get ε_X = 1.880, ε_Y = 4.573 and μ_Z = 0.612, and for α = 45° and kd = 0.25, we get ε_X = 1.848, ε_Y = 5.014 and μ_Z = 0.605. Also, in both cases, the value of the angle α is correctly retrieved (values are found to be respectively 0 and 44.62°).

Fig. 3 (a) Propagation of an incident plane wave (with θ = 45°) through a slab with kL = 5, containing a layered structure with kd = 0.5 and α = 0 (the layers are represented only on half of the slab to make visible the wavefield inside the slab). (b) Propagation through the slab (kL = 5) filled with the effective anisotropic medium defined by the retrieved parameters (α ≃ 0, ε_X = 1.880, ε_Y = 4.573 and μ_Z = 0.612).

Download Full Size | PDF

Fig. 4 Same representation as in Fig. 3 for (a) slanted layers in a slab with length kL = 2.5, with kd = 0.25 and α = 45°. (b) The effective parameters derived from the inversion method are α = 44.62°, ε_X = 1.848, ε_Y = 5.014 and μ_Z = 0.605.

Download Full Size | PDF

Comparisons are presented for α = 0 and kL = 5 (Fig. 3) and for α = 45° and kL = 2.5 (Fig. 4). As it can be seen, not only the scattering coefficients are correctly restituted (these scattering coefficients give the field outside the slab, and we find few percents of discrepancy on R and T in the presented cases), but also the fields inside the slab. The reader may notice however the scars visible in the wavefield obtained from direct numerics inside the layered structures (the layers are represented on only half of the slab to make visible the wavefield inside the slab).

To be more explicit, we now illustrate the different steps of our retrieval method in the case a = 45°. The preliminary step is illustrated in Fig. 5, for a ≠ 0, where the compensated transmission coefficient T_c is obtained from T (θ) and T (−θ), Eq. (21).

Fig. 5 Illustration of the preliminary step of the retrieval method: Calculation of the compensated transmission coefficient T_c, Eq. (21). Here T (θ) and T (−θ) have been calculated for slanted layers with α = 45° and kd = 0.5.

Download Full Size | PDF

Then, Figs. 6 and 7 show the retrieved index n, the impedance ratio ξ (step 1), and the parameter ε_y that is deduced from the previous two (step 2). In the presented cases, the sign ambiguity has been solved by selecting the sign of ξ̃ (from Eqs. (9) and (12), with cosθ > 0 and assuming ε_Y > 0) in order to impose the real part or the imaginary part of n to be positive; this corresponds to the selection of a right-going wave of the form e^inky with n either real positive for propagating waves (in the convention e^−iωt), or imaginary positive for an evanescent wave. Also, for our relatively low contrasts in mass density and bulk modulus, it has been sufficient to consider the principal determination of n in Eq. (12) (branch m = 0) and no branch ambiguity has appeared when varying the parameters θ, α and kd in the considered ranges.

Fig. 6 Step 1 of the retrieval method (main step): the impedance ratio ξ (θ) and the effective index (along the x-axis) n(θ) are retrieved, Eq. (12); for kd = 0.01 (○, blue), 0.1 (•, red) and 0.5 (⋄, green).

Download Full Size | PDF

Fig. 7 Step 2: ε_y is deduced from ξ and n; the final corresponding ε_X values are indicated. Same color and symbol conventions as in Fig. 6 are used.

Download Full Size | PDF

In Fig. 7, it can be observed that the retrieved results for ε_y calculated at different incidence angles agree well with each other for low kd with a small variation which increases with kd; this is expected since the assumption that the layered medium behaves as an effective medium fails when the wavelength becomes able to inspect the microstructure of size d. To anticipate, the values of ε_X (obtained after the whole retrieval procedure has been done) is indicated in each case in Fig. 7. As expected, it is angle and frequency independent, and it is found to be close to the homogenized one.

Finally, Fig. 8 illustrates the step 3 of the retrieval method; a fit of the quantity n²/ε_y is performed to determine μ_Z and the product ε_X ε_Y and Fig. 9 illustrates the step 4, where ε_xy is determined from T (θ)/T (−θ); here, a fit of T (θ)/T (−θ) has been performed in the form of a complex exponential function, according to Eq. (20).

Fig. 8 Step 3: Determination of μ_Z and of the product ε_X ε_Y by fitting the quantity n²/ε_y = μ_Z − sin² θ ε_y/ε_X ε_Y. Symbols correspond to the retrieved values (with the same color and symbol conventions as in Fig. 6), and plain lines to the fit.

Download Full Size | PDF

Fig. 9 Step 4: T (θ)/T (−θ) is fitted to find ε_xy, Eq. (20).

Download Full Size | PDF

The step 5 is straightforward: ε_x is deduced from Eq. (22), since ε_y and the product ε_X ε_Y have been retrieved. We find ε_x = 2.660, 2.664 and 2.750 for kd = 0.01, 0.1 and 0.5 respectively.

As previously said, alternatively to the components (ε_x, ε_y, ε_xy) of the permittivity tensor, we can calculate the full set of retrieved parameters (ε_X, ε_Y, μ_Z, α), Eqs. (23)–(24). They are given below

For kd = 0.01, ε_X = 1.845, ε_Y = 4.944, μ_Z = 0.607, α = 0.774,
For kd = 0.1, ε_X = 1.846, ε_Y = 4.953, μ_Z = 0.606, α = 0.775,
For kd = 0.5, ε_X = 1.855, ε_Y = 5.215, μ_Z = 0.603, α = 0.780.

Although we have not considered the extreme low frequency regime, the retrieved parameters are in quite good agreement with the homogenized parameters, Eq. (25), which are known to properly describe layered structures. This tends to validate our retrieval method. In the following section, we further inspect the retrieval method for varying α and kd.

4.2. Inspection of variations in α and k

Our retrieval method has been applied for varying α values, and kd = 0.01, 0.1 and 0.5. The results are shown in Fig. 10 for the set of parameters (ε_x, ε_y, ε_xy) and μ_Z. For α > 20°, a slab containing two unit cells ℓ_y = d/sinα has been considered. However, as α goes to zero, the length of the unit cell diverges; thus, we have imposed a fixed value of ℓ_y corresponding to α = 20°; this makes the results for α < 20° less reliable since the wave does not capture the periodicity of the structure along the y-axis, particularly for larger wavelength (this is why the results for kd = 0.5 and α < 20° have been omitted in the plots). For comparison, we have reported the behaviors of μ_Z and of (ε_x, ε_y, ε_xy) coming from homogenization theory, using the values given by Eq. (25) in Eq. (14). It can be seen that the retrieved parameters are close to the homogenized one, which is expected for these layered structures.

Fig. 10 Variation of the retrieved parameters with α, for kd = 0.01 (○, blue), kd = 0.1 (•, red) and kd = 0.5 (⋄, green). Gray dotted lines correspond to the homogenized values, Eqs. (25) and (14).

Download Full Size | PDF

Results for varying kd are reported in Fig. 11. To change the representations, we show (ε_X, ε_Y), the retrieved angle α_r, and μ_Z, and corresponding homogenized values are reported. Again, the retrieval values are close to the homogenized ones, with significant variations (sometime outside the plot) for kd approaching unity.

Fig. 11 Variation of the retrieved parameters (as defined in Eq. (6)) as a function of kd, for α = 0 (○, blue) and α = 45° (⋄, green). Gray dotted lines correspond to the homogenized values, Eqs. (25).

Download Full Size | PDF

5. Results for tilted elliptical particles and acoustic waves

As presented in Section 2, the case of acoustic waves is similar to the case of electromagnetic waves. In this section, we consider periodic arrangements of elliptical particles with relative (with respect to the surrounding medium) mass density ρ₁ = 100 and relative bulk modulus B₁ = 200. The particles are on a square lattice ℓ_x-periodic, and they are defined by the major and minor diameters a = 0.8ℓ_x and b = 0.2ℓ_x; being the case, they can be tilted by an angle α inside the unit cell (Figures 12).

Fig. 12 Square lattices of tilted elliptical particles; the unit cell is square ℓ_x × ℓ_x and (a, b) are the major and minor diameters. Note that tilting the ellipse inside the unit cell produces an array (b) with arrangement different from tilting the whole array (a). The surrounding medium has normalized mass density and bulk modulus; the elliptical particles have relative mass density ρ₁ and relative bulk modulus B₁.

Download Full Size | PDF

As first attempts, we consider the array of the Figure 12(b) with α = 45°. The retrieval method is performed starting from (R, T) calculated for ℓ_y = 4ℓ_x and kℓ_x = 10⁻²; we find α_r = 44.4° (close to value of the tilt angle α = 45°) and

ρ_{X} = 1.171, ρ_{Y} = 1.785, B = 1.138 .

These values are consistent with the results available in the literature for this configuration. Indeed, B is close to B^th, the volume average of the bulk moduli in the unit cell; with $φ = π a b / 4 ℓ_{x}^{2}$ the filling fraction of the particles, we have

B^{th} = {[φ + B_{1} + 1 - φ]}^{- 1} = 1.143,

often refereed as Wood’s law, see, e.g., [37]. For the effective density tensor, a simple expressions of (ρ_X, ρ_Y) is given in the dilute limit (small φ-values, and φ ≃ 0.126 in our case) [38]. This derivation is performed for electromagnetic waves and for an array corresponding to α = 0 (Eq. (38) in this reference); we give below the acoustic equivalent (accounting for the non convenient equivalences ε_X → ρ_Y and ε_Y → ρ_X):

{\begin{cases} ρ_{X}^{th} = \frac{2 - r_{X} - φ + ρ_{1} (r_{X} + φ)}{2 - r_{X} + φ + ρ_{1} (r_{X} - φ)} ≃ 1.168, \\ ρ_{Y}^{th} = \frac{2 - r_{Y} - φ + ρ_{1} (r_{Y} + φ)}{2 - r_{Y} + φ + ρ_{1} (r_{Y} - φ)} ≃ 1.853, \end{cases}

where we have defined r_X ≡ 2a/(a + b) and r_Y ≡ 2b/(a + b). The agreement with our retrieval values is reasonable (with differences less than 4%).

To get insight of the robustness of the retrieval parameters, we use the values of (ρ_X, ρ_Y, B) (in Eq. (26)) to characterize the effective behavior of arrays at significant higher frequency kℓ_x = 0.25, for a larger slab L = 50ℓ_x and for three α values (α = 0, 45°, 90°). Figures 13 show the fields calculated for an incident wave propagating through the real lattices, compared with the fields calculated for an incident wave propagating through the slabs filled with the effective anisotropic media; we used the effective parameters previously determined, Eq. (26) and the direction of anisotropy is fixed by the tilt angle α. The agreement is good in the three cases.

Fig. 13 Top panel: Propagation of an incident wave (with θ = 45°) through a slab with length L = 50ℓ_x containing tilted elliptical particles in a square lattice with unit cell of size kℓ_x = 0.25 (the ellipse is given by a = 0.8ℓ_x, b = 0.2ℓ_x). The particles are materialized with black line only on the lower part of the slab; the insets show a single particle in the unit cell. Bottom : Propagation through the slab of same length filled with the effective anisotropic medium defined using the retrieved parameters (ρ_X = 1.171, ρ_Y = 1.785, B = 1.138), and (a) α = 0, (b) α = 45°, α = 90°.

Download Full Size | PDF

The observed agreement, although intuitive, is not so trivial : indeed, the particles have been tilted by α inside the unit cell, but the lattice remains a square lattice; this is not the case if the whole array (at α = 0) is tilted of α; for instance, rotating the whole array of angle α = 45° would result in a hexagonal lattice while it remains a square lattice in our case. This suggests that the composition (ρ₁, B₁), the filling fraction φ and the orientations α of the particles are the pertinent parameters which characterize the effective medium, at least at leading order (as in the theoretical expressions Eqs. (27) and (28)).

To be more quantitative, we report in Figure 14 the variations of the retrieved parameters (ρ_X, ρ_Y, B, α_r) as a function of α; the retrieval method has been applied for kℓ_x = 10⁻² and with ℓ_y = Nℓ_x, N = 5, 15 and 30.

Fig. 14 Variation of the retrieved parameters with α for N = 5 (○, blue), N = 15 (•, red) and N = 30 (⋄, green). Dotted gray lines are the theoretical predictions, Eqs. (27) and (28).

Download Full Size | PDF

As a first comment, the retrieved parameters (ρ_X, ρ_Y, B) are found to be rather constant and the retrieved angle of anisotropic α_r is found to coincide with the tilt angle α within 2% for any N. This is consistent with the conclusions drawn from the fields in Figure 13, that the effective medium is not very sensitive to the type of lattice. Note also that (ρ_X, ρ_Y) satisfy the symmetry α → (180° − α), as expected.

Then, two facts are remarkable. (1) The amplitude in the variations of (ρ_X, ρ_Y) does not decrease with N, and the amplitude is quite significant for ρ_X (10 %). We also checked that it does not decrease for lower frequencies. We think that this variation may be attributable to the change in the lattice configurations when α varies, as previously commented. (2) An influence of N, the size of the array chosen for the inversion, has visible consequence. On the one hand, the retrieved bulk modulus B slightly increases with N. On the other hand, for ρ_X and ρ_Y, the symmetry α → (90° − α) is only approached for increasing N values. These two facts illustrate the “size dependance” of the effective anisotropic medium, which is a second order effect, compared with the dominant order (influence of (ρ₁, B₁), φ). Indeed the symmetry α → (90° − α) is expected for infinite lattices; for small N, it turns out that the effective parameters of an array with elongated particles orientated along x (square lattice with α = 0) or orientated along y (square lattice with α = 90°) are not identical. This is illustrated further in the last Figure 15 where we report the variations of the effective parameters with the array size N for α = 0 and α = 90°. While α_r is found independent of N (within 2% variation), the retrieved parameters ρ_X (and ρ_Y) for α = 0 and 90° tend to a commun value only asymptotically for large N.

Fig. 15 Variation of the retrieved parameters with the array size N (ℓ_y = Nℓ_x) for α = 0 (○, blue) and α = 90° (⋄, green). The dotted gray lines indicate the theoretical predictions, Eqs. (27) and (28).).

Download Full Size | PDF

6. Concluding remarks

We have proposed an extension of the well-known parameter retrieval method toward oblique incidence and arbitrary anisotropy directions of the effective medium. The retrieval method has been validated in a simple case of slanted layer structures, where the homogenization theory of layered medium is known to produce accurate results in the low frequency regime. Consistent results have been found for the retrieved parameters. We also reported results for tilted elliptical particles in the acoustic case; again, the retrieved parameters were consistent with available results of the literature, with dominant effects being due to the composition, the filling fraction and the orientation of the particles; beyond these dominant effects, the weak influence of the size of the medium has been evidenced.

However, in both cases, the effective parameters have been shown to be roughly independent on the characteristics of the incident wave (providing that the assumption of low frequency regime is reasonably satisfied). In many cases of interest, this is not the case; complex resonant structures are at least dependent on the frequency and may even depend on the incidence angle, and, in the worth cases, on the sample shape and size [23]. In these cases, the concept of effective medium becomes questionable or meaningless; enlightening critical discussions on this problem can be found in [21, 24, 25]. It is stressed in [25] that the retrieved parameters have to be understood rather as wave parameters able to restitute the scattering properties of the structure than as material parameters, independently of the scattering problem. Obvisouly, the interest in determining these effective parameters is to avoid heavy computations of a real macrostructure including its microstructure. Note that this implies at least the wave parameters to be independent on the size of the structures. With respect to this goal, it is of interest to introduce a conceptual homogeneous medium, associated to a known physics (its response with respect to wave propagation) and because the structures are involved, it is of interest that this medium possesses many degrees of freedom. This is the reason for previous improvements of the retrieval method and it is the meaning that we attribute to the extra degree of freedom considered in the present study (the angle α of the principal direction of anisotropy). Adding more and more degrees of freedom (within a prescribed known wave physics) should help to isolate more material parameters and less wave parameters. At least, this is what we can expect from such improved models.

Acknowledgments

The authors acknowledge the financial support of the Agence Nationale de la Recherche through the grant ANR ProCoMedia, project ANR-10-INTB-0914. A. C. and J-F. M. acknowledge the funding from the European Community’s Seventh Framework Program (FP7/2007–2013) under grant agreement no 285549: SIMPOSIUM project.

References and links

1. N. Engheta and R. W. Ziolkowski, eds. Metamaterials: Physics and Engineering Explorations (John Wiley & Sons, 2006). [CrossRef]

2. L. Solymar and E. Shamonina, Waves in Metamaterials (Oxford University, Oxford, 2009).

3. R. V. Craster and S. Guenneau, eds. Acoustic metamaterials: Negative Refraction, Imaging, Lensing and Cloaking, Vol. 166 (Springer, 2012).

4. J. Sánchez-Dehesa, D. Torrent, and J. Carbonell, “Anisotropic metamaterials as sensing devices in acoustics and electromagnetism,” Proc. of SPIE 8346, 834606 (2012). [CrossRef]

5. C.P. Berraquero, A. Maurel, P. Petitjeans, and V. Pagneux, “Experimental realization of a water-wave metamaterial shifter,” Phys. Rev. E 88, 051002(R) (2013). [CrossRef]

6. O. A. Oleinik, A. S. Shamaev, and G. A. Yosifian, Mathematical Problems in Elasticity and Homogenization (North Holland, Amsterdam, 1992).

7. T. Antonakakis, R. V. Craster, and S. Guenneau, “High-frequency homogenization of zero-frequency stop band photonic and phononic crystals,” New J. Phys. 15, 103014 (2013). [CrossRef]

8. P. Markos and C. M. Soukoulis, “Transmission properties and effective electromagnetic parameters of double negative metamaterials,” Opt. Express 11, 649–661 (2003). [CrossRef] [PubMed]

9. D. R. Smith, D. C. Vier, T. Koschny, and C. M. Soukoulis, “Electromagnetic parameter retrieval from inhomogeneous metamaterials,” Phys. Rev. E 71, 036617 (2005). [CrossRef]

10. X. Chen, T. M. Grzegorczyk, B. Wu, J. Pacheco, and J. A. Kong, “Robust method to retrieve the constitutive effective parameters of metamaterials,” Phys. Rev. E., 70, 016608 (2004). [CrossRef]

11. V. Fokin, M. Ambati, C. Sun, and X. Zhang, “Method for retrieving effective properties of locally resonant acoustic metamaterials,” Phys. Rev. B 76, 144302 (2007). [CrossRef]

12. D. R. Smith, S. Schultz, P. Markos, and C. M. Soukoulis, “Determination of effective permittivity and permeability of metamaterials from reflection and transmission coefficients,” Phys. Rev. B 65, 195104 (2002). [CrossRef]

13. H. Chen, J. Zhang, Y. Bai, Y. Luo, L. Ran, Q. Jiang, and J. A. Kong, “Experimental retrieval of the effective parameters of metamaterials based on a waveguide method,” Opt. Express 14, 12944–12949 (2006). [CrossRef] [PubMed]

14. Z. Liang and J. Li, “Extreme acoustic metamaterial by coiling up space,” Phys. Rev. Lett. 108, 114301 (2012). [CrossRef] [PubMed]

15. Y. Xie, B.-I. Popa, L. Zigoneanu, and S. A. Cummer, “Measurement of a broadband negative index with space-coiling acoustic metamaterials,” Phys. Rev. Lett. 110, 175501 (2013). [CrossRef] [PubMed]

16. B.-I. Popa and S. A. Cummer, “Design and characterization of broadband acoustic composite metamaterials,” Phys. Rev. B 80, 174303 (2009). [CrossRef]

17. L. Zigoneanu, B.-I. Popa, and S. A. Cummer, “Design and measurements of a broadband two-dimensional acoustic lens,” Phys. Rev. B 84, 024305 (2011). [CrossRef]

18. L. Zigoneanu, B.-I. Popa, A.F. Starr, and S. A. Cummer, “Design and measurements of a broadband two-dimensional acoustic metamaterial with anisotropic effective mass density,” J. Appl. Phys. 109, 054906 (2011). [CrossRef]

19. Y. Hao and R. Mittra, FDTD Modeling of Metamaterials: Theory and Applications, Chap. 7, (Artech House, Boston, 2009).

20. C. Menzel, T. Paul, C. Rockstuhl, T. Pertsch, S. Tretyakov, and F. Lederer, “Validity of effective material parameters for optical fishnet metamaterials,” Phys. Rev B, 81, 035320 (2010). [CrossRef]

21. S. B. Raghunathan and N. V. Budko, “Effective permittivity of finite inhomogeneous objects,” Phys. Rev. B 81, 054206 (2010). [CrossRef]

22. C. R. Simovski and S. A. Tretyakov, “On effective electromagnetic parameters of artificial nanostructured magnetic materials,” Photonics and Nanostructures - Fundamentals and Applications 8, 254–263 (2010). [CrossRef]

23. A. I. Căbuz, A. Nicolet, F. Zolla, D. Felbacq, and G. Bouchitté, “Homogenization of nonlocal wire metamaterial via a renormalization approach,” J. Opt. Soc. Am. B, 28, 1275–1282 (2011). [CrossRef]

24. C. R. Simovski, “On electromagnetic characterization and homogenization of nanostructured metamaterials,” J. Opt. 13, 013001 (2011). [CrossRef]

25. C. Menzel, C. Rockstuhl, T. Paul, F. Lederer, and T. Pertsch, “Retrieving effective parameters for metamaterials at oblique incidence,” Phys. Rev. B 77, 195328 (2008). [CrossRef]

26. F.-J. Hsieh and W.-C. Wang, “Full extraction methods to retrieve effective refractive index and parameters of a bianisotropic metamaterial based on material dispersion models,” J. Appl. Phys. 112, 064907 (2012) [CrossRef]

27. Z. H. Jiang, J. A. Bossard, X. Wang, and D. H. Werner, “Synthesizing metamaterials with angularly independent effective medium properties based on an anisotropic parameter retrieval technique coupled with a genetic algorithm,” J. Appl. Phys. 109, 013515 (2011) [CrossRef]

28. S. Arslanagić, T. V. Hansen, N. A. Mortensen, A. H. Gregersen, O. Sigmund, R. W. Ziolkowski, and O. Breinbjerg, “A review of the scattering-parameter extraction method with clarification of ambiguity issues in relation to metamaterial homogenization,” IEEE Antennas and Propagation Magazine 55, 91–106 (2013). [CrossRef]

29. T.-C. Yang, Y.-H. Yang, and T.-J. Yen, “An anisotropic negative refractive index medium operated at multiple-angle incidences,” Opt. Express 17, 24189–24197 (2009). [CrossRef]

30. B. L. Wu, W. J. Wang, J. Pacheco, X. Chen, J. Lu, T. M. Grzegorczyk, J. A. Kong, P. Kao, P. A. Theophelakes, and M. J. Hogan, “Anisotropic metamaterials as antenna substrate to enhance directivity,” Microw. and Opt. Techn. Lett. 48, 680–683 (2006). [CrossRef]

31. Y. Yuan, L. Shen, L. Ran, T. Jiang, J. Huangfu, and J. A. Kong, “Directive emission based on anisotropic meta-materials,” Phys. Rev. A 77, 053821 (2008). [CrossRef]

32. T. Jiang, Y. Luo, Z. Wang, L. Peng, J. Huangfu, W. Cui, W. Ma, H. Chen, and L. Ran, “Rainbow-like radiation from an omni-directional source placed in a uniaxial metamaterial slab,” Opt. Express 17, 7068–7073 (2009). [CrossRef] [PubMed]

33. A. Maurel, S. Félix, and J.-F. Mercier, “Enhanced transmission through gratings: Structural and geometrical effects,” Phys. Rev. B 88, 115416 (2013). [CrossRef]

34. T. M. Grzegorczyk, Z. M. Thomas, and J. A. Kong, “Inversion of critical angle and Brewster angle in anisotropic left-handed metamaterials,” Appl. Phys. Lett. 86, 251909 (2005). [CrossRef]

35. A Maurel, J.-F. Mercier, and S. Félix, “Wave propagation through penetrable scatterers in a waveguide and though a penetrable grating,” J. Acoust. Soc. Am. 135, 165–174 (2014). [CrossRef] [PubMed]

36. E. Popov and M. Nevière, “Grating theory: New equations in Fourier space leading to fast converging results for TM polarization,” J. Opt. Soc. Am. A 17, 1773–1784 (2000). [CrossRef]

37. D. Torrent and J. Sánchez-Dehesa, “Effective parameters of clusters of cylinders embedded in a nonviscous fluid or gas,” Phys. Rev. B 74, 224305 (2006). [CrossRef]

38. N. A. Nicorovici and R. C. McPhedran, “Transport properties of arrays of elliptical cylinders,” Phys. Rev. E 54(2), 1945–1957 (1996). [CrossRef]

Electromagnetics (TM)	Electromagnetics (TE)	Acoustics

Magnetic field H = (0, 0, H)	Electric field E = (0, 0, E)	pressure field p

If α = 0, 1/ε_Y	1/μ_Y	1/ρ_X
If α = 90°, 1/ε_X	1/μ_X	1/ρ_Y
Otherwise (arbitrary α), $\frac{1}{ε_{x}} = \frac{{sin}^{2} α}{ε_{X}} + \frac{{cos}^{2} α}{ε_{Y}}$	$\frac{1}{μ_{x}} = \frac{{sin}^{2} α}{μ_{X}} + \frac{{cos}^{2} α}{μ_{Y}}$	$\frac{1}{ρ_{x}} = \frac{{cos}^{2} α}{ρ_{X}} + \frac{{sin}^{2} α}{ρ_{Y}}$
$\frac{1}{ε_{y}} = \frac{{cos}^{2} α}{ε_{X}} + \frac{{sin}^{2} α}{ε_{Y}}$	$\frac{1}{μ_{y}} = \frac{{cos}^{2} α}{μ_{X}} + \frac{{sin}^{2} α}{μ_{Y}}$	$\frac{1}{ρ_{y}} = \frac{{sin}^{2} α}{ρ_{X}} + \frac{{cos}^{2} α}{ρ_{Y}}$
$\frac{1}{ε_{x y}} = - \frac{sin 2 α}{2} (\frac{1}{ε_{X}} - \frac{1}{ε_{Y}})$	$\frac{1}{μ_{x y}} = - \frac{sin 2 α}{2} (\frac{1}{μ_{X}} - \frac{1}{μ_{Y}})$	$\frac{1}{ρ_{x y}} = - \frac{sin 2 α}{2} (\frac{1}{ρ_{X}} - \frac{1}{ρ_{Y}})$

μ_Z	ε_Z	1/B

Generalized method for retrieving effective parameters of anisotropic metamaterials

Abstract

1. Introduction

2. 2D Wave equation in anisotropic media for electromagnetic and acoustic waves

3. Retrieval method

3.1. Case of anisotropy direction α = 0

3.2. Case of tilted principal axes, α ≠ 0

4. Results for layered dielectric media

4.1. Validation of the retrieval method

4.2. Inspection of variations in α and k

5. Results for tilted elliptical particles and acoustic waves

6. Concluding remarks

Acknowledgments

References and links

Cited By

Figures (15)

Tables (1)

Equations (28)

Optics Express