## Abstract

Over the past several years, metamaterials have been introduced and rapidly been adopted as a means of achieving unique electromagnetic material response. In metamaterials, artificially structured—often periodically positioned—inclusions replace the atoms and molecules of conventional materials. The scale of these inclusions is smaller than that of the electromagnetic wavelength of interest, so that a homogenized description applies. We present a homogenization technique in which macroscopic fields are determined via averaging the local fields obtained from a full-wave electromagnetic simulation or analytical calculation. The field-averaging method can be applied to homogenize any periodic structure with unit cells having inclusions of arbitrary geometry and material. By analyzing the dispersion diagrams and retrieved parameters found by field averaging, we review the properties of several basic metamaterial structures.

© 2006 Optical Society of America

## 1. INTRODUCTION

There has been a surge of interest over the past several years in the development of artificially structured periodic materials whose electromagnetic properties can be succinctly expressed in terms of homogenized material parameters. Unlike photonic crystal structures, whose typically complex electromagnetic modes are best described using the language of band structure and similar concepts associated with photonic crystals, the detailed description of an artificially structured composite can frequently be replaced by a small set of averaged parameters, such as the electric permittivity, magnetic permeability, refractive index, and wave impedance.

It might seem that this simplified description of an otherwise complex medium could result in all the interesting electromagnetic features of the medium being lost. On the contrary, work on bianisotropic and chiral materials over the past decade has shown the important role that artificial materials can play in terms of achieving unusual electromagnetic phenomena never before realized or difficult to obtain in conventional materials.[[1], [2], [3]] So successful have been the recent advances in artificial materials research that a new term has emerged and been adopted by the community to describe these novel structures: metamaterials.[[4]]

What constitutes a metamaterial is still not well defined, or at least the term is not uniformly applied by all authors; other well-established terms such as effective media and artificial dielectrics—even photonic crystals—could be considered to be referring to the same structures.[[5]] As conceived, the term metamaterial carried a somewhat different implication relative to previous terms for artificial structures: the derived electromagnetic response of a metamaterial is unachievable or unavailable in conventional materials.

Although the term metamaterial remains nebulous, the artificial structures that have become associated with the metamaterials typically have several features in common: metamaterials are described with homogenized material parameters; they are often based on conducting, resonant inclusions; the inclusions are positioned periodically; and the scale of periodicity is smaller—usually by a factor of 5 to 10 times—than that of the free-space wavelength of excitation. Structures that bear these features do not fit easily into the existing categories of artificial materials. For example, although the unit-cell dimension of the recently introduced metamaterials is less than the wavelength, it is only marginally so; the more relevant parameter ${k}_{0}d$—the free-space wavenumber multiplied by the unit-cell dimension—is of the order of unity. Homogenization theories are typically valid when the unit-cell size is insignificant with respect to the wavelength (the zero-frequency limit) and thus might be expected to result in a poor description of metamaterials.

Artificial structures possessing negative index and all the features noted above have become perhaps the best known examples of metamaterials. The possibility of a material with negative refractive index was hypothesized in 1968 by Veselago,[[6]] who showed that a negative refractive index would be a key property of a hypothetical composite whose electric permittivity and magnetic permeability were both less than zero. Although such a material condition is, in principle, achievable in a naturally occurring material, the fundamental atomic or molecular processes that lead to the electromagnetic response of materials do not appear to favor the overlapping electric and magnetic resonances that would be necessary to achieve negative index.[[7]]

In 2000, a metamaterial, based on conducting wires and split-ring resonators (SRRs), was demonstrated to have a negative refractive index over a region of microwave frequencies.[[8]] Given the absence of any conventional material with this property, the negative-index composite is consistent with the original and most strict definition of a metamaterial.[[4]] Since the introduction of this structure, the conducting wire[[9]] and SRR[[10]] inclusions and their variations have become widely used as the fundamental building blocks in artificial materials to provide tailored artificial electric and magnetic responses.[[11], [12], [13], [14], [15], [16], [17]]

Conductors are especially useful in the construction of artificial materials, as they strongly scatter electromagnetic fields. Moreover, conductors have an inherently inductive response due to the currents that flow within them, a response absent in structures composed of dielectric materials only. Thus, in metallic and metallodielectric metamaterials, the existence of both inductive and capacitive elements within the unit cell can result in strong resonances that couple either electrically or magnetically to the electromagnetic field. Being related to the local charge and current distributions rather than the overall geometry, these resonances can occur in structures much smaller than the free-space wavelength. Conducting inclusions are thus ideal as the basis of artificial structures with unique homogenized electromagnetic properties.

Given the recent successes that metamaterials have achieved and the hope that metamaterials will lead to enhanced electromagnetic devices and new applications, there is a need for accurate and efficient methods of determining the electromagnetic parameters for a given structure. The goal of assigning electromagnetic parameters to a composite or mixture of two or more component materials has long been the aim of traditional analytic homogenization theories, which have resulted in many well-known approximate mixing formulas, such as the Clausius–Mossotti and Maxwell–Garnett formulas, or the coherent potential approximation.[[18]] All known analytical methods, however, are valid under certain limitations and particular geometries or classes of structures. For metamaterials comprising conducting and possibly resonant elements, and for which the periodicity is not necessarily negligible relative to the free-space wavelength, analytical homogenization techniques are unreliable or not applicable.

Given that Maxwell’s equations can now be solved rapidly and efficiently—at least for domains relatively small with respect to the wavelength—with custom or commercial software running on modern personal computers, the local electric and magnetic fields, current distributions, and other quantities of interest can be found easily for virtually any type of structure. A purely numerical approach to the homogenization of artificial structures is thus feasible and avoids the limitations associated with prior analytic homogenization models.

One approach to the assignment of electromagnetic parameters to a simulated structure is via *S*-parameter retrieval.[[19]] In the *S*-parameter retrieval method, the complex reflection and transmission coefficients are calculated for a finite thickness of metamaterial, usually one unit cell, and compared with analytical formulas for the reflection and transmission coefficients of a homogeneous slab with the same thickness. The retrieval method can be applied to both simulated or measured *S*-parameter data, although thin samples (in the propagation direction) are best, and the retrieved electromagnetic parameters are not unique.

Despite its versatility and ease of use, *S*-parameter retrieval does not provide significant physical insight into the nature of the artificial material. The *S*-parameter observables—reflection and transmission coefficients—are not easily connected to analytical homogenization theories, which instead usually involve analysis of the local-field and current distributions.

In deriving the effective-medium properties of several periodic artificially magnetic structures, Pendry *et al.* presented a straightforward method of homogenization based on averaging the local fields according to[[10]]

**E**and

**H**fields, and $\overline{\mathbf{D}}$ and $\overline{\mathbf{B}}$ are obtained from integrals of the local

**D**and

**B**fields over the surfaces of the unit cell. The length of a side of the unit cell, assumed cubic for this example, is

*d*.

Equations (1) imply that a discrete set of averaged fields are to be defined at each unit cell, replacing the local fields that vary throughout the unit cell. In this sense, the averaging scheme of Eqs. (1) resembles a finite differencing of Maxwell’s equations, with a resulting spatial discretization analogous to those used in finite-difference time-domain and similar algorithms.[[20], [21]] In the case that $d\u2aa1\lambda $, where *λ* is the wavelength in the effective medium defined by the periodic structure, then we expect Eqs. (1) to provide an appropriate macroscopic description, from which we can ascertain the effective electromagnetic parameters.

The condition that *d*—effectively the scale of the inhomogeneity—be negligible with respect to the wavelength is equivalent to the statement that we are in the electrostatic or magnetostatic limit. In the static limit, one can formulate Eqs. (1) directly by considering the electrostatic or magnetostatic potentials.[[22]] In either case, the validity of Eqs. (1) would appear to be limited to the static regime. However, for most of the metamaterials demonstrated to date, *d* is not negligible with respect to the wavelength; yet, in these cases, the description of the structure as a homogeneous medium has produced sensible and useful results, implying that it should be possible to modify or extend the homogenization concepts embodied in Eqs. (1). In 2000, such a modification was suggested and shown to produce self-consistent results for the permittivity and permeability of wire and SRR media.[[23]]

Our intent here is to justify a numerically based homogenization scheme based on Eqs. (1), in which the local fields computed for one unit cell of a periodic structure are averaged to yield a set of macroscopic fields. Once having computed the macroscopic fields, we can then determine the constitutive relationships between the macroscopic fields, arriving at the effective electromagnetic parameters. We will see that there will be virtually no restrictions on the contents of the unit cell, nor will the unit cell necessarily need to be small in comparison with the wavelength.

## 2. FIELD-AVERAGING METHOD

To facilitate our development of a generalized metamaterial homogenization theory, we must first solve Maxwell’s equations subject to periodic boundary conditions to find the local fields associated with each propagating mode. We will restrict our discussion to the transverse modes, so that we need consider only the Maxwell curl equations:

Here, $\mathbf{E},\mathbf{B},\mathbf{D}$, and**H**are all local fields. Locally,

**D**and

**H**are related to

**E**and

**B**by the appropriate constitutive relationship for a particular homogeneous material. We assume the unit cell contains two or more such homogeneous materials such that the unit cell is inhomogeneous. In general, for a periodic medium for which

**a**is a primitive lattice vector, the solutions to Eqs. (2) are Bloch functions having the form

**q**is the Bloch or reciprocal-lattice vector.

As can be seen from Eq. (3), the boundary condition on the fields has the form

**q**that lie outside the range $-\pi \u2a7d\mathbf{q}\bullet \mathbf{a}\u2a7d+\pi $ can always be related to a

**q**value within this range having the same boundary condition. Frequency solutions for all

**q**values can thus be folded back into this first Brillouin zone, with the resulting band structure providing a specification of the electromagnetic mode structure of the computed unit cell.

We can compute the exact frequencies and local fields of the electromagnetic modes of a periodic structure by first specifying a value of **q** (or $\mathbf{q}\bullet \mathbf{a}$) and solving Eqs. (2) subject to the boundary conditions of Eq. (4). This method of solution is known as the eigenfrequency method, as the wave equation assumes the form of an eigenvalue equation once a value of **q** has been selected. For each value of **q** the eigensolver will return an infinite set of modes, each of which is continuous as a function of **q**.

Note from the form of the solutions in Eqs. (3, 4) that if we are uninterested in the details of the fields within a unit cell, then the Bloch-wave solutions throughout the periodic medium are similar to plane-wave solutions if we restrict our evaluation of the fields to a set of grid points that lie, for example, at the vertices or edges of the unit cells. By averaging Maxwell’s equations for the local fields over a unit cell, we arrive at a modified version of Maxwell’s equations in which the averaged fields are defined only at specific grid points. Thus, the exact electromagnetic response and geometric details of the inclusions within a unit cell are eliminated, encapsulated in the newly defined averaged (or macroscopic) fields $\overline{\mathbf{E}},\overline{\mathbf{B}},\overline{\mathbf{D}}$, and $\overline{\mathbf{H}}$. If we assume these macroscopic fields satisfy a constitutive relationship, then electromagnetic parameters such as $\overline{\epsilon}$ and $\overline{\mu}$ can be defined. We expect the homogenized material parameters to be consistent with other properties of the medium, such as the dispersion relation, $\mid \mathbf{q}\mid ={n}_{\mathrm{eff}}\omega \u2215c=\sqrt{\overline{\epsilon \mu}}\omega \u2215c$.

To illustrate the method by which we define field averages, consider the (unshaded) cube shown in Fig. 1 , each side of which has a length $2d$. The origin of the coordinate system is assumed to be at the center of the sphere, which represents the contents of one repeated unit cell of the infinite structure. Within the unit cell there is no restriction on the materials or geometry of the inclusions, so that we can generally allow for such properties as finite or infinite conductivity (perfect electric conductor) as well as arbitrary values of the local *ε* and *μ*. Within the cube, and exterior to the inclusions, the local ${E}_{x},{E}_{y}$, and ${B}_{z}$ fields are related by the first of the Maxwell curl equations in Eqs. (2):

Integrating the first of Eqs. (5) over the $x=+d$ surface, the second over the $y=+d$ surface, and the third over the $z=+d$ surface of the unit cell and applying Stokes’s theorem to the left-hand sides, we obtain

Using the definitions in Eqs. (7), we can replace the components of the curl equation, Eqs. (6), with

To solve for the electromagnetic fields in a structure, we must obtain a second equation for the averaged fields. Following the path that led to Eqs. (8), we start with the second Maxwell equation relating the curl of the magnetic field to the electric field, or

Writing down the components of Eq. (9) explicitly, we have**H**to those of

**D**as follows:

Equations (8, 13) form a set of finite-difference equations that can be solved to yield plane-wave solutions. To illustrate the nature of these solutions, consider a plane wave propagating along the *y* axis, polarized such that the electric field is directed along the *z* axis while the magnetic field is directed along the *x* axis. All fields are expected to have the spatial dependence $\mathrm{exp}(i\mathbf{q}\bullet \mathbf{x})$, where **x** is the position vector. From Eqs. (8) we find

It is instructive to consider the explicit form for ${\overline{\epsilon}}_{x}$ and ${\overline{\mu}}_{z}$ when the unit cell is empty. Using the constitutive relation between $\overline{\mathbf{E}}$ and $\overline{\mathbf{D}}$, for example, we relate the macroscopic fields to ${\overline{\epsilon}}_{x}$ and ${\overline{\mu}}_{z}$ as follows:

Equations (19, 20) reveal that the material parameters found from the field-averaging method exhibit spatial dispersion; that is, the material parameters are functions of the propagation vector. Inserting Eqs. (19, 20) into Eq. (17), we recover the correct dispersion relation for free space, or

In short, those terms dependent on the wave vector in the effective material parameters are exactly canceled by similar terms in the finite-difference dispersion relation of Eq. (17).Spatially dispersive medium parameters are less intuitive and less convenient to apply. However, the terms dependent on **q** in Eqs. (19, 20) are seen to be a result of the finite differencing or discretization of Maxwell’s equations and have no other physical origin. Since we recover the correct free-space dispersion relation, Eq. (21), via this scheme there is the hint that the averaged material parameters may have validity if we simply remove the term $\mathrm{sin}\left({q}_{y}d\right)\u2215{q}_{y}d$. When the unit cell is not empty, this procedure will obviously not lead to an exact result but may still be a reasonable approximation.

Consider the most general form for the local fields of the periodic system, as in Eq. (3). Assuming the same polarization and direction of propagation as above, the electric field for the periodic system has the form

where $u\left(y\right)$ is a periodic function that satisfies $u(y+2d)=u\left(y\right)$. Using Eq. (22) in Eq. (18), we have*y*. In a similar manner,

## 3. NUMERICAL EXAMPLES

#### 3A. Field Averaging Using Eigensolutions

Any method by which expressions or values for the local fields are obtained can be used in conjunction with the field averages [Eqs. (7, 12), for example] to obtain the effective electromagnetic material parameters. To illustrate the method here, we make use of the eigensolver in HFSS (Ansoft), a commercial finite-element-based electromagnetic simulation software package, to solve for the local fields in a unit cell of arbitrary composition. In the eigensolver module, periodic boundary conditions are applied to the bounding faces of the unit cell to simulate an infinite structure. The surfaces along one of the axes are related by periodic boundary conditions with a phase advance that varies between 0° and 180°. At each phase advance, a set of eigenvalues (mode frequencies) and eigenmodes (field distributions) are generated.

When the unit cell possesses a mirror symmetry, electric or magnetic boundary conditions can be substituted for the periodic boundary conditions. For example, if we seek the propagating modes corresponding to a unit cell with mirror symmetries about the planes perpendicular to the propagation direction, then applying electric boundaries along one of the axes and magnetic boundaries along the other will correctly simulate the effective medium, having the added benefit that the polarization of the solutions will be selected.

The advantage of the eigensolver method is that, by one’s examining the modes as a function of phase advance, the electromagnetic properties of the medium can be understood in the context of its photonic band structure. The disadvantage of the eigensolver is that it does not easily provide information in the bandgap regions—where *ε* or *μ* (but not both) are negative, for example, or where material losses are large. These difficulties can be avoided by one’s applying other numerical computational approaches, such as frequency- or time-domain simulators.

To obtain the field averages from the simulation data, we define integration lines along the edges of the cell to compute $\overline{\mathbf{E}}$ and $\overline{\mathbf{H}}$ and integration surfaces over the faces of the cell to compute $\overline{\mathbf{D}}$ and $\overline{\mathbf{B}}$. The integrals are performed using the macro capabilities of HFSS in the postprocessor. For the flux integrals ($\overline{\mathbf{D}}$ and $\overline{\mathbf{B}}$) taken over surfaces whose normals are perpendicular to the propagation direction, a factor of $\mathrm{sin}\left(qd\right)\u2215qd$ is removed from the average. For structures with no magnetoelectric coupling, we can relate the effective electromagnetic material parameters directly to the appropriate ratios of the averaged fields.

In the absence of losses, the fields and field ratios should be real; however, depending on the convention of the software used, the line integrals may contain a factor of $\mathrm{exp}(i\mathbf{q}\bullet {\mathbf{x}}_{0}-i\omega t)$, where ${\mathbf{x}}_{0}$ is the position in the cell along the propagation direction where the integration lines are defined and *t* is the time. To accommodate for this possibility, it is necessary to perform the field averages at two phases 90° apart, which can arbitrarily be represented as the real and imaginary parts of the fields. In the absence of loss in the system and with no magnetoelectric coupling, the 0° and 90° field averages should maintain a constant ratio as a function of phase advance, such that multiplying the complex fields by a phase factor will result in purely real fields.

In certain geometries, such as a medium of continuous wires, a conducting element may pass continuously from one unit cell to the next. For the flux averages to be correct in this case, the formulation of the effective-medium parameters given above would need to be modified so as to include the contribution of the current flowing between the cells. As an alternative, a thin slice could be made in the conductor at the surface so that the current is no longer continuous; if the slice is very small, the electromagnetic properties of the simulated structure will not be affected. The field averaging would then include the depolarizing fields within the gap, effectively accounting for the contribution of the current in the conductor.

To avoid having to account for elements that transition between unit cells when we compute the effective material parameters, we can often apply an alternative approach that avoids the flux averages. Because the eigensolver provides mode frequencies as a function of the phase advance *ϕ*, or wavenumber $k=\varphi \u2215d$, it is straightforward to determine the effective index *n* as the ratio of a given eigenfrequency to the wavenumber, or $n=\omega \u2215k$. Once having determined the index, we can characterize the material by computing the wave impedance $z=\overline{E}\u2215\overline{H}$. In the absence of magnetoelectric coupling the relations

#### 3B. Dielectric Sphere Medium

Before considering metallic metamaterial structures, it is useful to first consider a system that can also be analyzed using analytical effective-medium theories. Such a system is an array of dielectric spheres, as shown in Fig. 2 (left). When the spheres are small relative to their spacing, the Maxwell–Garnett equation should yield the effective dielectric of the composite, or[[18]]

The computed permittivity and permeability corresponding to four different sphere radii are shown in Fig. 3 as a function of $d\u2215\lambda $ and reveal several interesting features. At zero $d\u2215\lambda $, the recovered permittivity is in exact agreement with that predicted by Eq. (30), shown as the horizontal dashed lines. The permeability approaches unity at zero frequency, as would be expected in the electrostatic limit. At finite wavelengths, however, the computed values of the permittivity and the permeability deviate characteristically from the static limit, with the composite permeability taking values less than unity and the composite permittivity taking values greater than the static value. That the material parameters disperse as a function of frequency should not be considered surprising, since a bandgap is approached (by definition) when the phase advance approaches 180°. At the band edge the index approaches a fixed value set by the lattice constant, $n=\mathrm{c}\pi \u2215\left(d\omega \right)$, but the permittivity and permeability can have resonant forms. (Note that we now define the unit-cell length as *d* rather than $2d$.) This phenomenon has been reported in studies of the effective-medium parameters of metamaterials[[25]] and will be seen to play an even greater role in the metamaterials to be discussed.

#### 3C. Wire Medium

As a first example of a conducting metamaterial, we consider a one-dimensional array of wires with square cross section, as illustrated in the inset to Fig. 4 . When the electric field is polarized along the wire and the wave propagates in the plane perpendicular to the wire axis, then the wire medium can be reasonably described—at least at lower frequencies—by an effective permittivity having the form

The effective plasma frequency, ${\omega}_{p}$, is a constant related to the lattice parameters, such as the unit-cell length of the wire medium and the cross-sectional area of the wire. Losses, ignored here, since the simulations assume perfect conductors, add a damping term in Eq. (31).Equation (31) was reported in 1996[[9]] and has since been confirmed theoretically by numerous other authors.[[26], [27], [28], [29], [30]] The striking feature shown by Eq. (31) is that the effective permittivity is negative for frequencies below ${f}_{p}={\omega}_{p}\u2215\left(2\pi \right)$. An eigenmode simulation of the wire medium produces the dispersion curve plotted in Fig. 4, which reveals a single propagating band that starts at ${f}_{p}$, in general qualitative agreement with Eq. (31).

The simple plasmonic form of Eq. (31) excludes effects due to the periodicity of the medium. In particular, the propagating band shown in Fig. 4 terminates at a second band of forbidden frequencies at the band edge (where $k=\pi \u2215d$). Because the effects due to the higher-frequency band structure are not included in Eq. (31), we expect the properties of the actual wire-medium structure to diverge from the simple model as the higher frequencies are approached. The material parameters as determined via the field-averaging technique are presented in Fig. 5 (solid curves). For comparison, the permittivity for the homogeneous plasmonic medium [Eq. (31)] is plotted as the dashed curve and can be seen to deviate significantly from the wire medium at frequencies above ${f}_{p}$. The presence of the upper-frequency bandgap forces the permittivity of the wire structure toward zero, with the permeability approaching a resonance. Clearly, since the wavenumber is fixed at the band edge to be $k=\pi \u2215d$, the effective index must also be fixed as $n=ck\u2215\omega $; thus, if one of the material components approaches a resonance, the other must approach zero such that the square root of their product is equal to the effective index at the band edge.[[25]]

The dispersion curves exhibited by the actual wire structure, as shown in Fig. 5, are nonintuitive, so it is reasonable to question whether the effective material parameters produced are in any way meaningful. The effective material parameters should, for example, describe the manner in which a plane wave is reflected and transmitted through a finite section of the composite. We can assess this by comparing the material parameters retrieved via field averaging with those retrieved by an *S*-parameter inversion. Using an algorithm described previously to invert the transmission and reflection coefficients obtained from an *S*-parameter simulation on the same wire structure,[[19]] we obtained the results for the permittivity and permeability shown as the gray or black dots in Fig. 5. As the figure reveals, there is excellent agreement between the field averaging and the *S*-parameter retrieval methods, with any discrepancy within the numerical error for these calculations. Although we do not show similar comparisons in the subsequent figures, we find generally that the material parameters retrieved by the two methods consistently agree to within the numerical accuracy of our simulations.

The results of Figs. 4, 5 suggest a relatively simple picture of the homogenized properties of the anisotropic wire medium. In particular, Eq. (31) is a valid description of the medium so long as incident waves are restricted to propagate perpendicular to the axes of the wires, as modeled in the simulations presented here. For a wave whose propagation vector makes an arbitrary angle relative to the wire axes and that has a component of its electric field along the wire axes, the resulting mode structure is considerably more complex. For such a wave, Eq. (31) no longer holds; the transverse permittivity displays strong spatial dispersion, and additional longitudinal modes may also be excited.[[26], [27], [28], [29], [30]] Although we do not discuss the analysis of such modes here, we note that the boundary conditions can be adjusted in the eigensolver to simulate propagation along nonprincipal axes of the wire medium, with the field-averaging method subsequently applied.

#### 3D. Split-Ring Resonator Medium

The wire medium provides a conceptually straightforward means of achieving a basic Drude or Drude–Lorentz type of electric response. Because of this, the wire medium in various incarnations has achieved prominence in metamaterial designs. In the same way, the conducting split-ring-resonator (SRR) medium introduced in 1999 by Pendry *et al.* has become the archetypical magnetic component in metamaterials.[[10]] Applying an effective-medium theory based on Eqs. (1), Pendry *et al.* showed that the electromagnetic response of the SRR medium–an example of which is shown in Fig. 6a –could be approximately described by the resonant form

*F*is a constant relating to the filling factor of the unit cell. Because the SRR is composed of conductors with no inherently magnetic components, the magnetic response of the SRR can be considered artificial, the result of induced circulating currents.

To illustrate the properties of the SRR medium, we consider the structure shown in Fig. 6a and in the insets to Figs. 7, 8, 9. To facilitate the discussion of the SRR properties, we assume a set of axes fixed to the SRRs as shown in Fig. 6a. Analytical calculations have revealed that the electromagnetic material parameters of the SRR medium have the approximate form

with all other diagonal elements equal to unity, all other off-diagonal elements equal to zero, and the constants $\alpha ,\beta ,\gamma $, and*δ*related to the particular SRR geometry.[[31]]

Although introduced as a means of obtaining a magnetic response, Eqs. (33) show that the SRR is actually a highly anisotropic, electromagnetically complex structure, exhibiting not only a resonant magnetic response but also resonant electric and magnetoelectric responses.[[32]] The resonant magnetoelectric coupling, in particular, results from the lack of a mirror plane symmetry about the axis lying in the SRR plane and perpendicular to the ring gaps. A medium arranged as shown in Fig. 6b thus exhibits bianisotropy. One can form medium with the magnetoelectric coupling term equal to zero by alternating the SRRs so that gaps alternately face each other and away from each other, as shown in Fig. 6c. To remove the resonant electric response as well, one must alter the inherent symmetry of the SRR; a structure that minimizes both the resonant electric and the magnetoelectric couplings is proposed in Fig. 6d.

Consider first the SRR medium of Fig. 6b, oriented as shown in the leftmost inset to Fig. 7 . The wave is assumed to propagate along the *y* axis, with the electric field polarized along the *z* axis. Under this circumstance the macroscopic fields are related by

The mode frequency as a function of phase advance is shown as the black curve in Fig. 7, which reveals to good approximation the resonant form predicted by assuming that the SRR medium of Fig. 6 is approximated by a permeability like that of Eq. (32). The lower branch rises from zero frequency to the band-edge frequency $\left({\omega}_{m0}\right)$, after which a frequency gap occurs, followed by a second propagation band. The computed dispersion curve deviates from that predicted by Eq. (32) owing to the presence of the second frequency gap. The retrieved permeability, shown in Fig. 8 as the black curve, reveals the expected resonant form: unity at zero frequency; large positive values near the resonant frequency; no propagating modes found in the gap region (where the permeability is negative); and zero at the start of the next propagating band and increasing with frequency toward the next gap.

The effective-medium arguments leading to Eqs. (32, 33) also predict that the SRR medium, for the same polarization of incident wave, will have a nontrivial electric response. Pendry *et al.* argued that an approximate form for the ${\epsilon}_{zz}$, valid at low frequencies, should be $1\u2215\sqrt{1-F}$, with $F\u2a7d1$.[[10]] This correction is necessary, as it maintains the light line at or below *c*. The permittivity, plotted as the gray curve in Fig. 8, is indeed approximately constant at low frequencies with magnitude significantly larger than unity. Somewhat surprising, however, is that the permittivity approaches zero at the resonance and appears to have an antiresonant form in general. As described above, the antiresonant nature of the permittivity in this case has been discussed as a consequence of the underlying periodicity of the medium.[[25]]

Figure 8 reveals the expected result that the SRR medium can respond with a resonant permeability and, in particular, can provide usable negative values of the permeability from which negative-index materials can be formed. However, this interpretation implies a particular relationship between the direction and the polarization of the exciting wave and the SRR. In Fig. 7, a second dispersion curve is shown, for which the exciting fields are such that a predominantly electric response is expected. In this case, the relevant constitutive relationship is

*x*axis in this case) must vanish, and only the values of ${\epsilon}_{yy}$ and ${\mu}_{zz}$ are needed to characterize the SRR medium.

Figure 9 shows the retrieved material parameters for the SRR oriented as in the inset to Fig. 9. It is found that ${\epsilon}_{yy}$ exhibits a strong resonance, as expected. Although the analytic theory leading Eqs. (33) predicts no response for ${\mu}_{zz}$, we see there is actually an antiresonance in ${\mu}_{zz}$, again owing to the finite size of the unit cell.

#### 3E. Calculation of the Bianisotropy Parameters

In the two orientations of the SRR relative to the incident wave considered in Subsection 3D, consideration of the magnetoelectric tensor elements could be avoided. There may be structures designed to exhibit magnetoelectric coupling, however, where it is desired to calculate the magnitude of the coupling. If we orient the SRR as shown in the inset to Fig. 10 , the constitutive relations have the form

For a lossless system, the material parameters and the magnetoelectric coupling terms are expected to be real-valued functions. An inspection of Eqs. (38) thus shows that the field ratios, denoted by $\zeta ,\xi $, and *z*, will now be complex-valued functions. The various products and ratios in Eqs. (38, 39) must therefore be performed using the complex values obtained from the 0° and 90° field averages, as described earlier. In general, we find that stability is much more difficult to achieve for structures that possess magnetoelectric coupling, especially near the resonant frequency of the composite structure.

In Fig. 10 we compare the significant material parameters of the SRR shown in Fig. 6a for frequencies up to the bandgap. Equations (38) were used to calculate ${\kappa}_{yx}$, with $\gamma ,\xi $, and *z* being calculated from field averages with the SRR in the bianisotropic configuration and ${\epsilon}_{yy}$ and ${\mu}_{xx}$ taken from the field averages obtained for the SRR in the electric and magnetic configurations, respectively. The behavior of all terms is consistent with analytic predictions.[[31], [32]]

#### 3F. Negative-Index Medium

As described above, the wire medium possesses a frequency region over which $\epsilon <0$ for waves propagating perpendicular to the wire axes and polarized such that their electric field vectors lie in the direction of the wires. Likewise, the SRRs possess a frequency region over which $\mu <0$, for the SRR medium properly oriented as described above. If the two structures are combined such that the SRR-medium magnetic resonance overlaps with the $\epsilon <0$ frequency region of the wire medium, then the composite medium should exhibit negative index. The structure is similar in concept to the negative-index composite first reported in 2000.[[8]]

The inset to Fig. 11 shows a combined wire-and-SRR unit cell. The wire has a square cross section, with side length 0.04 of the unit-cell length. The SRR and unit cell have the same dimensions as described above in Subsection 3D. The lower computed dispersion curve in Fig. 11 shows the familiar signature of negative index, with the phase and group velocities having opposite sign.

Applying the field-averaging procedure on the eigenfields corresponding to the passbands, we recover the effective-medium parameters shown in Fig. 12 . For modes in the lowest passband, the relevant components of *ε* and *μ* are seen to be negative, as is the recovered index. Also, as expected, the recovered material parameters for the upper band are all positive.

## 4. CONCLUSION

Field averaging provides a useful means of characterizing the effective-medium parameters of an arbitrary metamaterial. It is particularly useful for complex structures in which the transmitted and reflected amplitudes are not easily interpreted; in such cases it can be difficult to choose the appropriate solution (out of the available set of solutions) that matches the expected effective-medium theory. In contrast, the material parameters obtained by the field-averaging method are not subject to the multiple branches that occur in the *S*-parameter method. When combined with dispersion curve calculations, field averaging provides a wealth of information distinct and complementary to that obtained by *S*-parameter analysis and retrieval.

The calculations presented here were performed using a script written in the Visual Basic Scripting language that runs within the HFSS software. The script steps through phase advances between 0° and 180°, computing the eigenmodes and eigenfields and performing the field averages at each phase step. Steps in intervals of 20° were used in all the simulations presented above. Input parameters are provided in an EXCEL spreadsheet read by the script, with the computed mode frequencies and field averages being output to the same EXCEL spreadsheet. In this way, the process of characterizing a metamaterial can be automated, allowing also for sophisticated iterative design cycles.

The examples presented are typical of metamaterials, in that the size of the unit cell is generally smaller but by no means negligible compared with the wavelength, as indicated by the normalized frequency scale in all the figures. Still, the field-averaging method provides physically reasonable retrieved parameters, which are also consistent with those obtained by *S*-parameter retrieval. The results of both methods indicate that metamaterials can be usefully—even if only approximately—described by effective-medium parameters. It can be readily seen from the examples presented that somewhat less intuitive material parameters are obtained as the higher-order band structure (above the first Brillouin zone) is approached. Although the field-averaging method appears to provide consistent results at these frequencies, it can be expected to be unreliable at frequencies where multiple passbands are present.

Finally, we should note that we have presented examples for which the unit cell is relatively simple, consisting of only one or two elements, generally designed to have desired properties in only one or two spatial dimensions. When the unit cell becomes more complicated, having multiple elements, additional bands are found in the eigenmode simulation that make the analysis more difficult. These bands may not be excited in a plane-wave transmission–reflection experiment but do introduce additional complexity into the eigenmode spectrum. For these more intricate structures, the field-averaging method can, in principle, be applied, but the analysis becomes more involved.

## ACKNOWLEDGMENTS

The authors acknowledge informative discussions with S. Tretyakov and R. Marqués. This work was supported by a Multidisciplinary University Research Initiative, sponsored by Defense Advanced Research Projects Agency, Office of Naval Research (grant N00014-01-1-0803).

Corresponding author D. R. Smith may be reached by e-mail at drsmith@ee.duke.edu.

**1. **M. M.I. Saadoun and N. Engheta, “A reciprocal phase shifter using novel pseudochiral or omega-medium,” Microwave Opt. Technol. Lett. **5**, 184–188 (1992). [CrossRef]

**2. **I. V. Lindell, A. H. Sihvola, S. A. Tretyakov, and A. J. Viitanen, *Electromagnetic Waves in Chiral and Bi-Isotropic Media* (Artech House, 1994).

**3. **W. S. Weiglhofer and A. Lakhtakia, eds., *Introduction to Complex Mediums for Optics and Electromagnetic* (SPIE, 2003). [CrossRef]

**4. **R. M. Walser, “Electromagnetic metamaterials,” in Complex Mediums II: Beyond Linear Isotropic Dielectrics, A. Lakhtakia, W. S. Weiglhofer, and I. J. Hodgkinson, eds., Proc. SPIE **4467**, 1–15 (2001). [CrossRef]

**5. **A. Sihvola, “Electromagnetic emergence in metamaterials,” in *Advances in Electromagnetics of Complex Media and Metamaterials*, S. Zouhdi, A. Sihvola, and M. Arsalane, eds., Vol. 89 of NATO Science Series II: Mathematics, Physics, and Chemistry (Kluwer Academic, 2003), pp. 3–17. [CrossRef]

**6. **V. G. Veselago, “The electrodynamics of substances with simultaneously negative values of *ε* and *μ*,” Sov. Phys. Usp. **10**, 509–514 (1968). [CrossRef]

**7. **J. B. Pendry and D. R. Smith, “Reversing light with negative refraction,” Phys. Today **57**, 37–43 (2004). [CrossRef]

**8. **D. R. Smith, W. Padilla, D. C. Vier, S. C. Nemat-Nasser, and S. Schultz, “A composite medium with simultaneously negative permeability and permittivity,” Phys. Rev. Lett. **84**, 4184–4187 (2000). [CrossRef] [PubMed]

**9. **J. B. Pendry, A. J. Holden, W. J. Stewart, and I. Youngs, “Extremely low frequency plasmons in metallic mesostructures,” Phys. Rev. Lett. **76**, 4773–4776 (1996). [CrossRef] [PubMed]

**10. **J. B. Pendry, A. J. Holden, D. J. Robbins, and W. J. Stewart, “Magnetism from conductors and enhanced nonlinear phenomena,” IEEE Trans. Microwave Theory Tech. **47**, 2075–2084 (1999). [CrossRef]

**11. **P. Markoš and C. M. Soukoulis, “Numerical studies of left-handed materials and arrays of split ring resonators,” Phys. Rev. E **65**, 036622 (2002). [CrossRef]

**12. **M. Bayindir, K. Aydin, E. Ozbay, P. Markoš, and C. M. Soukoulis, “Transmission properties of composite metamaterials in free space,” Appl. Phys. Lett. **81**, 120–122 (2002). [CrossRef]

**13. **C. G. Parazzoli, R. B. Greegor, K. Li, B. E.C. Koltenbah, and M. Tanielian, “Experimental verification and simulation of negative index of refraction using Snell’s law,” Phys. Rev. Lett. **90**, 107401 (2003). [CrossRef] [PubMed]

**14. **P. F. Loschialpo, D. L. Smith, D. W. Forester, F. J. Rachford, and J. Schelleng, “Electromagnetic waves focused by a negative-index planar lens,” Phys. Rev. E **67**, 025602 (2003). [CrossRef]

**15. **A. A. Houck, J. B. Brock, and I. L. Chuang, “Experimental observations of a left-handed material that obeys Snell’s law,” Phys. Rev. Lett. **90**, 137401 (2003). [CrossRef] [PubMed]

**16. **R. W. Ziolkowski, “Design, fabrication, and testing of double negative metamaterials,” IEEE Trans. Antennas Propag. **51**, 1516–1529 (2003). [CrossRef]

**17. **F. Falcone, T. Lopetegi, M. A.G. Laso, J. D. Baena, J. Bonache, M. Beruete, R. Marques, F. Martin, and M. Sorolla, “Babinet principle applied to the design of metasurfaces and metamaterials,” Phys. Rev. Lett. **93**, 197401 (2004). [CrossRef] [PubMed]

**18. **V. M. Shalaev, “Electromagnetic properties of small-particle composites,” Phys. Rep. **272**, 61–137 (1996). [CrossRef]

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

**20. **K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. **14**, 302–307 (1966). [CrossRef]

**21. **J. B. Pendry, “Calculating photonic band structure,” J. Phys. Condens. Matter **8**, 1085–1108 (1996). [CrossRef]

**22. **D. J. Bergman, “The dielectric constant of a composite material—a problem in classical physics,” Phys. Rep. **9**, 377–407 (1978). [CrossRef]

**23. **D. R. Smith, D. C. Vier, N. Kroll, and S. Schultz, “Direct calculation of permeability and permittivity for a left-handed metamaterial,” Appl. Phys. Lett. **77**, 2246–2248 (2000). [CrossRef]

**24. **J. A. Kong, *Electromagnetic Wave Theory*, 2nd ed. (Wiley, 1990), Sec. 2.6.

**25. **T. Koschny, P. Markoš, D. R. Smith, and C. M. Soukoulis, “Resonant and antiresonant frequency dependence of the effective parameters of metamaterials,” Phys. Rev. E **68**, 065602 (2003). [CrossRef]

**26. **S. I. Maslovski, S. A. Tretyakov, and P. A. Belov, “Wire media with negative effective permittivity: a quasi-static model,” Microwave Opt. Technol. Lett. **35**, 47–51 (2002). [CrossRef]

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

**28. **A. L. Pokrovsky, “Analytical and numerical studies of wire-mesh metallic photonic crystals,” Phys. Rev. B **69**, 195108 (2004). [CrossRef]

**29. **C. R. Simovski and P. A. Belov, “Low-frequency spatial dispersion in wire media,” Phys. Rev. E **70**, 046616 (2004). [CrossRef]

**30. **M. G. Silveirinha and C. A. Fernandes, “Homogenization of 3-D-connected and nonconnected wire metamaterials,” IEEE Trans. Microwave Theory Tech. **53**, 1418–1430 (2005). [CrossRef]

**31. **R. Marqués, F. Medina, and R. Rafii-El-Idrissi, “Role of bianisotropy in negative permeability and left-handed metamaterials,” Phys. Rev. B **65**, 144440 (2002). [CrossRef]

**32. **R. Marqués, F. Mesa, J. Martel, and F. Medina, “Comparative analysis of edge- and broadside–coupled split ring resonators for metamaterial design—theory and experiments,” IEEE Trans. Antennas Propag. **51**, 2572–2581 (2003). [CrossRef]