Micro- and macroscopic simulation of periodic metamaterials

In order to characterize three-dimensional, lefthanded metamaterials (LHM) we use electromagnetic field simulations of unit cells. For waves traveling in one of the main directions of the periodic LHM-arrays, the analysis is concentrated on the calculation of global quantities of the unit cells, such as scattering parameters or dispersion diagrams, and a careful interpretation of the results. We show that the concept of equivalent material values – which may be negative in a narrow frequency range – can be validated by large “global” simulations of a wedge structure. We also discuss the limitations of this concept, since in some cases the macroscopic behavior of an LHM cannot be accurately described by equivalent material values.


Introduction
In his triggering publication in 1968, Veselago proposed the existence of left-handed waves for materials with double negative material parameters, i.e. a permittivity and permeability below zero (Veselago, 1968).During the last 15 years, a number of realizations of such materials have been developed as composites: periodic arrays of unit cells with metallic or dielectric substructures on a microscopic (sub-wavelength) scale.The reversal of many concepts of wave propagation, such as phase velocity or the angle of refraction, and the possibility to engineer the material parameters within a large and unusual range, open many prospectives for the design of new components in microwave techniques and optics.
Being composites with local fields on a microscopic scale, the negative material parameters of LHM realizations can only be understood in terms of effective material values, and their determination for a given setup has become a central is-Correspondence to: R. Schuhmann (schuhmann@tet.upb.de)sue of the research on this field.Clearly, these effective material values will show dispersion, and very often the wellknown Lorentz-and Drude formulas from atomic material models can be successfully applied.
The extraction methods used today either start from a microscopic point of view, i.e. local fields within the unit cells, or from global quantities such as reflection and transmission coefficients of single cells or of larger arrays.The raw data can be obtained by analytical considerations (Pendry et al., 1999) or measurements (Shelby et al., 2001), or -with increasing importance -by field simulations (Smith et al., 2000;Weiland et al., 2001).Typical results are the complex effective permittivity and permeability values as functions of frequency.
The effective material values which are obtained by such techniques are an important result for the synthesis of the composites.However, their main purpose is to be used in the design process of large-scale components, where the novel properties of metamaterials shall be profitably used.Thus, the extraction results should of course be as accurate as possible, in the sense that the equivalent model can describe the behavior of the metamaterial within the interesting frequency band.However, the model should also be physical, i.e. it should be passive (as long as no active materials are present on the microscopic scale) and last but not least it should be causal.As a critical review of published results (such as, e.g., in Simovski and Tretyakov, 2007) shows, these properties are sometimes not fulfilled, and it seems as if this is not a result of shortcomings in the respective implementations.
The rest of the paper is organized as follows: first, we shortly introduce the simulation model for a typical metamaterial unit cell, which is used as an example to discuss various extraction methods for effective material values.In Sect.2.4 we describe a validation approach for this extraction process and show first results of this validation for the so-called PFDM approach.Finally, in Sect. 3 we propose to characterize metamaterials by the more general (and of Published by Copernicus Publications on behalf of the URSI Landesausschuss in der Bundesrepublik Deutschland e.V. R. Schuhmann et al.: Micro-and macroscopic simulation of periodic metamaterials course well-known) approach of dispersion diagrams in order to avoid the problems with non-physical effective material parameters.Using such dispersion diagrams as reference solution it is shown that the characterization of metamaterials based on scattering matrices can be enhanced, if higher order modes at the ports are considered.

Numerical model
We concentrate on 3-D metamaterial concepts such as the well-known split ring resonator (SRR) and similar structures, which can be analyzed by electromagnetic field simulation.All simulations are performed using commercial as well as laboratory implementations of the Finite Integration Technique, FIT (Weiland, 1977(Weiland, , 1996)), a full-wave simulations method based on a three-dimensional computational grid.If used on Cartesian grids and in time domain, in combination with the leapfrog integration scheme, the FIT is computationally equivalent to the FDTD method (Yee, 1966).However, the more general formulation also allows for various other solvers, e.g. for the eigenvalue formulation in frequency domain with periodic boundary conditions, as discussed below.The implementations include the described extraction methods for effective materials as well as all excitation techniques and boundary conditions needed for the simulation models.Further, for the so-called macroscopic simulations advanced material models are available to treat dispersive permittivity and permeability parameters in time domain (Gutschling et al., 2000).The structure modeling and parts of the postprocessing and visualization are performed by the commercial tool CST Microwave Studio (MWS, by Computer Simulation Technology GmbH, 2007).

Extraction methods
One of the first ideas to define effective material values in metamaterials goes back to Pendry et al. (1999).The local field distribution in three-dimensional unit cells of SRRs and similar structures is estimated by quasi-static approximations, and the effective material values are defined as the ratio of averaged field quantities.Following this technique, the dispersion of the effective permeability is found to be of Lorentz type.For the effective permittivity of an array of conducting wires a similar investigation leads to a Drudetype dispersion curve, and both together may constitute a double-negative material.
In order to consider also more general geometric shapes of the microscopic elements and their cross-coupling effects, several implementations of the field averaging approach based on simulated field data have been published, e.g. in Smith et al. (2000); Weiland et al. (2001).However, this type of methods have turned out to be not very robust, e.g. the desired dispersion types can not always be confirmed.
The second (and most often used) class of methods are oriented on classical measurement techniques and try to define the effective medium based on scattering matrices (SM).Typically these are the reflection and transmission coefficients of plane waves illuminating a metamaterial array with finite length.The effective material parameters can be determined by inverting the dependency of typical quantities (like S-parameters S 11 , S 21 , or the wave impedance Z) on and µ.This step is performed at each frequency to gain the desired (ω) and µ(ω)-curves, respectively.Again, the expected Lorentz-and Drude-dispersion models are often not revealed.
This last problem can be healed by the so-called parameter fitting of dispersive models (PFDM) approach in Lubkowski et al. (2006), where these dispersion models are assumed to hold a-priori, and their parameters are fitted to the measured or simulated SM data.Thus, the resulting material values are physical (i.e.passive and causal) by definition, whereas the achievable accuracy and the valid bandwidth are not easy to determine a-priori.The example given below shows, however, that the PFDM approach can be a very useful tool, which opens a straight-forward way for a validation of the effective LHM model.

Properties and limitations of equivalent material models
A main issue of effective material values is the question whether they are simultaneously negative in order to support left-handed waves.However, like all dispersive models, the permittivity and permeability are complex quantities, making the question of sign a bit more involved.A nice overview on this topic has recently been published in Tsukerman (2007).
The next issue to be discussed is the passivity and causality of the effective material, representing important constraints in the calculation of the complex permittivity and permeability over frequency.From above mentioned extraction methods, only the analytical (quasi-static) estimation and the PFDM method yield provably physical results.On the other side, the standard SM-based techniques typically face a couple of problems: First, they often suffer from numerical relicts within the extraction process such as numerical noise in the raw data, ill-conditioned inversion formulas, or improperly chosen branches of complex roots.Next, such methods only consider one principle wave traveling through the metamaterial (in most cases a plane wave).This issue is discussed in Sect.3.2.
As a result, the material curves may show imaginary parts (i.e.losses) in lossless simulation models, the wrong sign of such imaginary parts (i.e. an active effective material), or a violation of the Kramers-Kronig relations and thus of the principle of causality.
Finally, all scattering matrix based results (including those from PFDM) are only valid for one single direction of propagation of the waves.From the PFDM approach we obtain exact Drude and Lorentz models, respectively.The double-negative frequency band is between 9.7 GHz and 10.2 GHz.

Validation of equivalent material models
In spite of the large number of publications on effective material values in metamaterials, not many of them include a thorough validation of the results.Concerning the extraction of the material values from local fields, a validation should compare the behavior of the full (microscopic) LHM model with the effective one.Such a setup is shown in Fig. 1 (full geometric data are given in Lubkowski et al., 2006Lubkowski et al., , 2007)): Following the idea of the famous wedge experiment published by Shelby et al. (2001), the propagation of a wave through a wedge-shaped LHM and its refraction at the interface to the surrounding air are studied.In a first simulation, all details of the LHM, i.e. the microscopic structure of an array of split ring resonators and wires are considered in a single, large numerical model.For the given wedge we need several millions of grid points (with six unknowns per point), and it is obvious that this procedure is not feasible for arbitrary components.The results are compared to a second simulation, where the LHM is represented by a homogeneous block of material with the dispersive material properties gained by the PFDM-extraction process.Note that this extraction is performed in advance, based on a separate (and computationally cheap) FIT simulation of one single unit cell of the array.The results from the PFDM extraction are the parameters of exact Drude-and Lorentz-type dispersion models for the permittivity and the permeability, respectively (cf.Fig. 2).They are fed into the advanced material models within FIT for the second wedge simulation.Thus, the effective model has the potential to save a considerable amount of computation resources, since its realization in the simulation tools is by far less expensive compared to the microscopic model.In our case, we are able to downscale the computational effort by approximately two orders of magnitude.
As an exemplary result, Fig. 3 shows the fields at one frequency in the left-handed frequency band (f =10 GHz): Obviously, a negative angle of refraction is supported in both simulations.Moreover, the setup also allows for a quantitative analysis of the wave propagation and thus of the extracted material parameters.A more detailed investigation will be presented in a subsequent publication.However, these preliminary results already clearly indicate that the PFDM extraction process can be sufficiently accurate, at least in some important parts of the desired frequency spectrum.It is also important to note, that for this type of validation the dispersion data of the effective materials are needed, e.g. as a direct result of the PFDM approach.In case of non-physical material values, violating the passivity or causality rules, the effective simulation will most probably become unstable.

Dispersion diagrams as reference solutions
One of the principle weaknesses of the standard SM-based approach is that it only considers plane waves both in the excitation of the unit cells and in the analysis of reflected waves.The accuracy of the material characterization can be improved, if also higher modes (in terms of the aperture of a unit cell) are taken into account.However, this idea hinders the definition of effective materials, since a simple homogeneous (equivalent) medium will not be able to represent all effects of a scattering matrix in higher dimensions.
In order to establish a reference solution for this type of analysis, the characterization of metamaterials using dispersion curves is used.Starting again with the field simulation of single unit cells of the arrays, such dispersion curves can be obtained in frequency domain using periodic boundary conditions with a varying phase angle ϕ (see also Smith et al., 2000).The corresponding phase factor exp(j ϕ) is imposed in the two longitudinal boundary conditions (the ports of the SM approach), and from the arising eigenvalue problem we obtain one point of each ω(ϕ)-curve per simulation run.

Higher order transfer matrices
An alternative way is to calculate the higher-order scattering matrices (S ij ), which relates the wave amplitudes a i , b i at the "ports" to each other, i.e. the normalized amplitudes of the plane wave and the higher order mode fields at the aperture of the unit cells.By imposing the Bloch condition (1) we can solve for the dispersion relation ϕ(ω).Note that for single modes at the ports (plane waves only) and the correlated 2×2-S-matrices, the expressions (with the vacuum material parameters 0 and µ 0 , the unit cell size in propagation direction d z and the characteristic impedance Z c ) convert this quantity back into the effective material parameters (Simovski and Tretyakov, 2007).
3.3 Quantitative impact of higher-order modes But we concentrate here on a more detailed analysis of the dispersion relation ϕ(ω).Figure 4 shows the dispersion diagram for two typical SRR-wire unit cells with different lengths (array parameters, cf.Fig. 5) and left-handed frequency bands around 4.8 GHz.Obviously there exists a discrepancy between the dispersion relation obtained by a fundamental mode scattering matrix and a 3-D full-wave eigenmode simulation with periodic boundary conditions in propagation direction, if the unit cell size d z is small.The field distribution of 3-D full-wave eigenmode simulation at the port plane in Fig. 5 reveals, that the Bloch wave existing in the lattice can not be modeled as a simple plane wave since the field in the aperture is contaminated by higher order modes.
To improve the agreement between the dispersion relations obtained by these two approaches, it is necessary to take into account higher order modes.The extraction of the desired ϕ(ω)-curve out of the simulated 2n×2n S-matrix (for n modes at each port) is performed by the intermediate step of the transfer scattering matrix formulation Note that all the internal quantities are matrices and vectors now.The formulas for the transformation of S-into Tmatrices can be found in standard textbooks.We now solve for eigenvalues of the T-matrix, which describe the desired Bloch waves.From the phase angles of the propagating ones (those which have a magnitude of 1), we arrive at an improved dispersion relation which includes now the higher order modes at the ports.A typical convergence behavior is shown in Fig. 6.As the number of considered modes increases, the deviation to the 3-D full-wave eigenmode simulation decreases.
As expected, the need for higher-order modes in the scattering matrix approach depends on the size of the metamaterial unit-cells.The fundamental Bloch wave can only be modeled as a plane wave at the port plane, if the ports are sufficiently far away from the SRR and wire elements, and the internal substructures are not too densely packed.Yet a densely packed substructure may be desirable, since the complete bandwidth of the LHM band is increased (cf.Fig. 4).
The analysis also shows, how the impact of higherorder modes can be quantitatively assessed by a comparison with the frequency-domain reference solution with periodic boundaries.In cases where no or only a small number of higher-order modes are necessary for accurate results, the standard time-domain SM approach (with broadband excitation signals) keeps its important advantage to be computationally more efficient than the multiple solution of eigenvalue problems.

Conclusions
The concept of effective material values for electromagnetic metamaterials is very popular, since there seems to be a simple relation between such parameters and the phenomenon of left-handed waves.We have shortly discussed the various extraction methods for effective materials, and we have shown, that in some cases their results can be validated by large scale field simulations.
However, many difficulties which have been reported in the extraction of effective material values from measured or simulated data seem to be not only due to insufficient implementations.The concept itself is rather limited in its practical application: First, such parameters are only valid for one specific polarization, one specific frequency and one specific angle of incidence.Second, as we have shown in our analysis, the situation gets even more complicated if the structure size gets too small.Some more limitations which we have not discussed here may appear if the dependence between electric and magnetic effects gets too strong, or if interfaces between such materials have to be properly modeled.
A more robust way to characterize such materials is the well-established technique of dispersion diagrams of periodic structures, where all such effects are inherently in-  cluded.Nevertheless, the concept of effective material parameter may be beneficial in some applications, if its limitations are kept in mind.

Fig. 1 .Fig. 2 .
Fig. 1.Validation of effective material models: Simulations setup of an LHM wedge with microscopic (left) and effective (right) simulation.The microscopic model contains all geometric details of a 16×9 array of SRR-type unit cells.

Fig. 3 .
Fig. 3. Results of the LHM-wedge simulations: Both models (microscopic, left, and macroscopic, right) show a negative angle of refraction and thus the desired left-handed behavior.

Fig. 5 .
Fig. 5. Electrical fields in port planes for three different lengths of the unit cell.(All two-dimensional field plots of three different structures are projected in one single figure.)

Fig. 6 .
Fig.6.Convergence behavior: Accuracy of the dispersion diagram depending on the number of considered modes in the higher-order scattering matrix approach.