1. Introduction

Metamaterials are by definition artificial structures composed of specifically tailored subwavelength unit cells that may collectively support exotic wave properties and interesting macroscopic phenomena, beyond what is available in any of their constituents at the frequencies of interest. Their development has advanced the possibility of material engineering by allowing a precise control of the wave interaction at the subwavelength scale. Metamaterial properties may overcome several limitations of conventional electromagnetic and optical devices, as in the case of antennas [1–4] and lenses [5–8].

It is challenging to predict the performance of a device that includes these artificial structures because of the complex scattering effects arising from a large number of subwavelength, and often resonant, inclusions forming arrayed structures. One approach to significantly reduce this complexity is based on homogenizing these arrays and defining effective parameters that can capture the bulk electromagnetic response of the structure. This approach is well established for mixtures and natural crystals, for which conventional techniques based on volume averaging have been developed over the years [9,10]. These mixing rules, however, are usually based on quasi-static approximations and ignore strong coupling effects among particles that are common in complex structures. Indeed, metamaterials offer challenges that are not found in simple collections of particles, natural materials, or mixtures: they are usually formed by densely packed arrays of resonant inclusions, whose strong collective interaction can produce spatial dispersion effects and other associated complications; moreover, the size and complex geometry of the inclusions do not allow neglecting the phase variation within each unit cell, which invalidates common quasi-static assumptions [11]. As a consequence, effective parameters retrieved for metamaterials often produce nonphysical results [12,13], such as violations of both energy conservation and Kramers-Kronig relations [14,15]. These issues recently inspired several researchers to reconsider the definition of metamaterial constitutive relations more rigorously [16–23].

Among these recent works, a homogenization model [22,23] has been put forward that defines macroscopic averaged quantities by taking into account the full dynamic and complex interaction within a dense periodic array of small inclusions near resonance. The main result of this work is a set of physically meaningful effective constitutive parameters that can accurately describe the metamaterial response as that of a bulk homogeneous material. It was shown in [22,23] that the homogenization model is independent of the form of excitation, which makes it appealing for applications in realistic device configurations. It was also shown that it becomes necessary to consider magneto-electric coupling effects even when simple center-symmetric inclusions with no bianisotropy or chirality are considered, since lattice effects can introduce inherent asymmetric coupling in the wave propagation, consistent with recent numerical findings [21]. This effect becomes especially important near the array resonances, where spatial dispersion can dominate the bulk response of metamaterials. When no sources are embedded in the metamaterial, the derived constitutive relations may be simplified by combining the magneto-electric parameters into a set of equivalent parameters that describe eigen-modal propagation along the array [23,24]. This leads to simplified constitutive relations without magneto-electric cross-coupling terms, but the associated equivalent parameters become inherently nonlocal even in the long-wavelength limit, may not satisfy the constraints dictated by passivity and causality conditions in local materials, and their dispersion features may produce wrong signs of their imaginary parts and negative slopes versus frequency. Nevertheless, these simplified parameters are consistent with the ones commonly retrieved using conventional techniques, like the Nicholson-Ross-Weir (NRW) method [12,13] and provide a compact and often adequate description of the scattering properties of the metamaterial.

The homogenization model in [22,23] was developed using infinitely periodic arrays of subwavelength inclusions and is generally valid for $k_{0}d<1$, where $k_0$ is the background wave number and d is the array period. In this paper, the validity of this model is investigated by analyzing scattering from finite metamaterial objects formed by three-dimensional (3-D) arrays of magnetodielectric spheres and comparing the scattered fields to those scattered from equivalent homogeneous objects that occupy essentially the same region in space as the array and the constitutive parameters of these objects are determined by using the homogenization model in [22,23]. When excited by identical sources, it is expected that the fields scattered by the metamaterial and the equivalent object will be only approximately equal because edge or truncation effects, which should diminish for larger arrays, are not taken into account in the homogenization model. In the following, we quantify the differences between the scattered fields and analyze their dependence on the artificial structure’s size and shape.

The objective of this paper is to determine the validity and accuracy of this homogenization model [22,23] for finite metamaterials. For this purpose, a numerical approach is adopted, where the metamaterial array size is gradually increased, the fields scattered by the array and its equivalent object are computed, and the results contrasted. For this validation approach to be useful, there are two critical requirements from the numerical method: (i) the scattered fields must be computed accurately in both the original and equivalent problem so that the differences between the simulation results can be attributed to the homogenization model and not to numerical errors, i.e., numerical errors must be much smaller than the errors due to inherent approximations introduced in the homogenization process; (ii) scattering from large 3-D arrays of subwavelength resonant inclusions and 3-D bulk objects with negative effective parameters must be computed. These requirements cause large computational costs, ill conditioning, and low-frequency breakdown for conventional numerical methods and general-purpose commercial software based on them. Thus, in this paper, a recently developed full-wave simulator that employs a well-conditioned surface-integral-equation formulation [25–27] and an FFT-accelerated iterative solution method is used [28–30]. Using this simulator, the near- and far-field response of increasingly larger 3-D cylindrical and cubic arrays of magnetodielectric spheres are compared to those of 3-D homogeneous objects with equivalent parameters determined by the homogenization model. The results show that, as the metamaterial size increases, the homogenized model can more accurately represent the array, and edge effects may indeed be neglected. The results also offer interesting physical insights into the effects of granularity at the edges of metamaterial arrays and into how the collective response of a periodic array can support negative phase velocity even in finite-size arrays of inclusions forming negative-index objects.

The rest of the paper is organized as follows: Section 2 briefly reviews the homogenization theory and the definition of effective parameters on which the modeling is based. Section 3 highlights the numerical challenges and reviews the methods applied to tackle them. Section 4 presents the application and validation of the homogenization method for a disk-shaped array in two different frequency regimes, with double-positive and double-negative effective properties, respectively. Section 5 focuses on a cube-shaped array with double-negative effective properties and investigates the effect of the array size on the accuracy of the homogenization method. Finally, we draw our conclusions in Section 6.

2. Homogenization theory for magnetodielectric metamaterial arrays

Consider a metamaterial composed of a 3-D array of spheres with radius $a$ arranged in a cubic lattice with period $d$. In the following, we assume that that the sphere permittivity is ${\epsilon }_d=13.8{\epsilon }_0$, permeability is ${\mu }_d=11.0{\mu }_0$, $d=25\text{mm}$, and $a=0.45d=11.25\text{mm}$. Two observations about these parameter choices are in order: (i) A relatively large ${\mu }_d$ is used for the inclusions in this paper even though the inherent magnetic response of materials in the optical regime is usually weak; this is mainly for the sake of consistency with previous papers, where these exact parameters were used to realize negative-index metamaterials at microwave frequencies [19–24]. The concepts presented in this paper may be directly applied to other optical metamaterial designs that realize magnetic effects with plasmonic inclusions [31,32], or also extended to purely electric effects. (ii) Material losses of the inclusions are neglected; this is in the interest of simplicity and to keep a clear boundary between propagation and bandgap region. We note that both the homogenization method and the numerical method used in this paper have been validated for inclusions with losses [23,29].

As shown in [24], the considered inclusions support electric and magnetic resonances at frequencies close to $k_{0}d\cong 0.74$, which implies that the array may strongly interact with the impinging wave in the long-wavelength regime. Under these conditions, the inclusions can be well represented by their dipolar response, which was used in [19] to predict a reasonable bandwidth over which negative-index propagation is expected along the array. It was further shown in [7] that in this frequency range the array possess quasi-isotropic properties and it may be a promising candidate to realize metamaterial lenses. After applying the homogenization theory [22], this array of center-symmetric, isotropic inclusions can be described with time-harmonic macroscopic constitutive relations (with suppressed $e^{i\omega t}$ dependence) of the form

In the absence of sources embedded in the array, the constitutive relations in Eq. (1) may be written in the more common form

Figure 1 shows the homogenized effective and equivalent parameters determined based on Eqs. (1) and (2) using the dipolar approximation to describe the inclusion interaction with the impinging wave [22]. The figure also shows the conventional NRW parameters [13,31] retrieved from a finite metamaterial slab consisting of eight unit cells and excited at normal incidence, after proper corrections to eliminate Fabry-Perot artifacts, as described in [33]. The number of unit cells, or the overall slab thickness, for the retrieval process was selected to be large enough so that end-effects from the slab boundaries may be neglected and small enough to allow a fast retrieval procedure. Far from the array resonance and the corresponding bandgap, marked by the yellow shade around $k_{0}d=0.67$, effective, equivalent, and NRW parameters nicely overlap with each other. Near the bandgap, however, equivalent and NRW parameters bend downwards towards unity, whereas the effective parameters grow following a typical Lorentzian line shape. A similar trend is found for the permittivity curves at the upper bandgap edge. Equivalent and corrected NRW parameters align very well, because there are no sources embedded in the metamaterial and therefore only the dominant eigenmode is excited [31, 33, 34]. This ensures that the dipolar approximation holds well in this frequency range, despite the higher-order multipolar coupling that may exist in a dense array. The effective magnetoelectric parameter ${\chi }_{eff}^o$ shown in Fig. 1(c) takes into account the nonlocal effects hidden in the equivalent parameters near the bandgap region, where they diverge from the effective ones. It is important to remark that this region is the one of most interest for metamaterial applications, due to the strong wave interaction that leads to negative effective and equivalent parameters, and a bulk negative index of refraction.

 

Fig. 1 The (a) permittivity, (b) permeability, and (c) magneto-electric coefficient that are used to homogenize an infinitely periodic metamaterial array composed of magnetodielectric spheres as specified in Section 2. Effective, equivalent, and NRW retrieval parameters are calculated using Eq. (1), Eq. (2), and the method in [31]. The dashed vertical lines in all plots indicate the frequencies of interest in the following numerical simulation.

Download Full Size | PDF

3. Numerical method

In order to investigate whether the homogenized parameters can successfully model finite metamaterial objects, a fast surface integral equation solver is used to simulate scattering from several sample finite arrays of the magnetodielectric spheres and the results are compared to those obtained from corresponding homogeneous objects with the equivalent parameters determined above. Classical surface integral equation solvers have appealing properties for scattering analysis [28]; however, the problems of interest are particularly challenging for three reasons: (i) large-scale arrays can contain thousands of magnetodielectric elements and millions of unknowns, (ii) homogenized objects have negative material parameters that have not been typically considered in conventional solvers, (iii) both scattering scenarios can lead to ill-conditioned matrices due to their resonant features. In this paper, a well-conditioned and FFT-accelerated surface integral equation solver is used to perform the necessary simulations. The solver employs combined field integral equations formulated for each piecewise homogeneous region to improve the convergence of iterative method-of-moments solution and the multiple-grid adaptive integral method to accelerate the computations at each iteration. Detailed formulation, validation, and complexity analysis of the solver can be found in [28–30]. It was shown in [30] that this solver can successfully address all the above challenges and can be reliably used to validate the homogenization model discussed in the previous section when applied to finite metamaterial objects. The following numerical simulations based on this approach analyze near-field distributions within and around the metamaterial structures and the far-field radar cross section (RCS). In all simulations, diagonal preconditioning was used and the iterative solutions were terminated when the L2 norm of the relative residual error was less than 10⁻³.

4. Near-field simulations of finite metamaterial cylindrical disks

We first consider a disk-shaped finite array composed of two parallel layers of the magnetodielectric spheres specified in Section 2, stacked along the $\widehat{x}$ axis. The spheres in each layer are arranged in a square lattice with period $d$, and the array is truncated to include only the spheres that are entirely contained within a circle of radius $12d$. The corresponding total number of spheres is 848 and the array is excited by an $\widehat{x}$-polarized plane wave propagating in the $-\widehat{z}$-direction at the frequencies $f_{\text{DPS}}$ (for which $k_{0}d=0.524$) and $f_{\text{DNG}}$ ($k_{0}d=0.761$), which correspond to double-positive and double-negative regimes, respectively, as indicated by the vertical dashed lines in Fig. 1. The equivalent homogeneous disk was chosen to be of radius $12d$ and height $1.9d$ to match the height of the sphere array. The computational costs for simulating the problem at both frequencies are compared in Table 1. Both geometries were meshed with an average edge length ~$1/10$ of the minimum wavelength, requiring $N=1404288$ and $N=10800$ edges for the array and the homogeneous disk, respectively. The total computational time is the sum of the matrix fill time and the product of the number of iterations with the solve time per iteration. It is obvious inspecting Table 1 that huge advantages are obtained in both computational time and memory usage by analyzing the homogenized model, as expected. It should be noted that the number of required iterations depends on the initial guess for the iterative solver as well as on the material parameters, and more advanced pre-conditioning techniques, not considered here, may significantly reduce this parameter. At frequency $f_{\text{DPS}}$, the required number of iterations for the metamaterial array is larger than for the homogeneous disk, whereas in the DNG regime at frequency $f_{\text{DNG}}$, the resonant DNG-DPS boundary on the rim of the homogeneous disk introduces a larger number of iterations [29]. Still, the total computational time remains much lower than for the array, due to the significantly reduced number of unknowns and solution time per iteration.

Table 1. Computational requirements for the simulations in Figs. 2-3.

View Table | View all tables in this article

At frequency $f_{\text{DPS}}$, the equivalent parameters of the object, extracted from Fig. 1, are ${\epsilon }_1=2.752{\epsilon }_0$ and ${\mu }_1=2.275{\mu }_0$. The corresponding electric field amplitude distribution in the central transverse plane between the two layers of spheres is plotted in Figs. 2(a) and 2(b), for the metamaterial array and the homogenized disk, respectively. The figure shows that the homogenized model qualitatively captures the complexity of the wave interaction in this regime, producing a similar near-field distribution around the object. Even the distribution inside the array is reproduced to a good degree by the homogeneous model in terms of spatial distribution and phase fronts, showing a focusing point inside the disk towards the rear end (indicated by black arrows in Fig. 2). Obviously, small gaps between neighboring inclusions produce larger field values in the array, which are not found in the homogenized model. Despite these differences, good qualitative agreement is found and the homogenized model appears to capture the physics of the problem in this long-wavelength scenario. Figure 2(c) presents the normalized difference of the field distributions for the metamaterial array and the homogeneous disk, relative to the impinging field. It can be seen that the most significant differences arise inside and at the back of the array. We should mention in this context that we should not expect that a homogenization model based on the dipolar approximation would be able to predict the complex near-field interaction between neighboring inclusions in a densely packed array. Better performance may be expected using the retrieved parameters obtained from the generalized NRW retrieval method based on this homogenization theory, as recently proposed in [35], but this goes beyond the scope of this work.

 

Fig. 2 Disk-shaped finite array in the double-positive regime. The disks are placed in the y-z plane. Electric near-field distribution (in dB scale) due to a plane wave traveling toward $-\widehat{z}$ for (a) the finite array and (b) the homogeneous disk with ${\epsilon }_1=2.752{\epsilon }_0$ and ${\mu }_1=2.275{\mu }_0$. (c) Normalized difference between the two fields. (d) Snapshot of the time-domain electric field as it propagates through the disk at $y=0$ for the disk-shaped array and homogeneous disk. See also Media 1 for time-domain animations of electric fields in the y-z plane.

Download Full Size | PDF

Figure 2(d) shows a time-domain snapshot of the electric field across the center of the disk [e.g., at $y=0$ in Figs. 2(a)–2(c)]. The field inside the metamaterial array experiences larger oscillation due to the strong field interaction between inclusions. This figure also shows a very consistent wave pattern entering the disk, larger differences inside the disk, and a similar wave pattern being eventually achieved in the shadow region after about three wavelengths. Animations in the supplemental material show the close agreement in the near-field time-domain electric field for the metamaterial array and the homogenized disk. Several minor ripples are seen in the metamaterial case; these can be attributed to the strong field interactions between finite-size spheres. The phase fronts in the shadow region do not perfectly match each other because of the small thickness of the metamaterial disk that brings significant edge effects to the field distribution. In the next section, we analyze a thicker metamaterial object that is shown to mitigate these edge effects and demonstrate better matching of wave fronts even at the back of the array.

The simulated frequency in Fig. 2 is far from the band with resonant effects; indeed, the equivalent bulk parameters are positive and essentially represent a volume-average permittivity and permeability. Nevertheless, the agreement is encouraging and shows the utility of the proposed homogenization approach. More interesting effects are expected in the resonant region. Figure 3 shows the results for a similar investigation in the double-negative regime at the frequency $f_{\text{DNG}}$. Homogenized permittivity and permeability in this case are both negative and correspond to ${\epsilon }_2=-0.9936{\epsilon }_0$ and ${\mu }_2=-0.6063{\mu }_0$, respectively. The field distribution obtained in this case for the homogeneous disk still captures the overall scattering response well, but it reproduces the details of the field distribution somewhat less accurately within the array, because $k_{0}d$ is larger and the inclusions are operating near resonance. The differences are also evident around the disk edge, where strong resonant fields are observed in Fig. 3(b). Figure 3(c) presents the difference of the field distributions for the metamaterial array and the homogeneous disk normalized to the impinging field. The main scattering effects are qualitatively well reproduced by the homogenized model and the errors are again mostly concentrated at the back of the array.

 

Fig. 3 Disk-shaped finite array in the double-negative regime. The disks are placed in the y-z plane. Electric near-field distribution (in dB scale) due to a plane wave traveling toward $-\widehat{z}$ for (a) the finite array and (b) the homogeneous disk with ${\epsilon }_2=-0.9936{\epsilon }_0$ and ${\mu }_2=-0.6063{\mu }_0$. (c) Normalized difference between the two fields. (d) Snapshot of the time-domain electric field as it propagates through the disk at $y=0$ for the disk-shaped array and homogeneous disk. See also Media 2 for time-domain animations of electric fields in the y-z plane.

Download Full Size | PDF

The comparison between near-field distributions is better visualized by the phase variation shown in Fig. 3(d). Again, like Fig. 2(d), the fields do not completely overlap, as expected, once it enters the resonant array of spheres. Nonetheless, the homogenized model predicts with good accuracy the field distribution right outside the disk. It is important to note that negative slope and backward phase propagation are not observed in these objects (see animations in supplemental materials) despite the double-negative effective parameters. This is because the metamaterial array and homogenized disk are too thin in the transverse dimension (the disk is composed of two layers of spheres) and edge effects dominate; in fact, additional simulations not presented here show that, as the thickness of the disk increases, backward phase propagation becomes more prominent in the array and is visible when the array is 10 spheres thick. Backward phase propagation in thicker and larger objects is explored further in the next section. Interestingly, the field distribution and the absence of backward-wave propagation are captured also in the homogenized model; the results agree well with the metamaterial array even for this scenario, which further validates this homogenization approach.

5. Near- and far-field simulations of metamaterial cubes

In this section we focus on cube-shaped metamaterial arrays, studying how the size affects the accuracy of the adopted homogenization model. The number of elements in the arrays is varied from 4 to 10 per side. This topology is selected because it is not subject to staircasing at the boundaries as in the previous geometry; this allows us to focus only on the size dependence avoiding the truncation issues associated with curved boundaries, which may be responsible for some of the minor discrepancies highlighted in the previous examples. For brevity we focus here only on the double-negative regime, which is most interesting and challenging in terms of proper homogenization and numerical simulation. The metamaterial cubes are formed by the same magnetodielectric sphere arrays as in the previous section and are excited by $\widehat{x}$-polarized plane waves propagating in the $-\widehat{z}$ direction at the same frequency $f_{\text{DNG}}$ as in Fig. 4, Fig. 5, Fig. 6.

 

Fig. 4 Far-field (a, b) and near-field (c-h) distributions at f_DNG = 1.454 GHz for a 4-element cube-shaped array and homogenized cube in the double-negative regime with ${\epsilon }_2=-0.9936{\epsilon }_0$ and ${\mu }_2=-0.6063{\mu }_0$. (i) Snapshot in time of the electric field component $E_x$ as it propagates toward the $-\widehat{z}$ direction through the center of the cubes ($y=0$). See also Media 3 for time-domain animations in the y-z plane.

Download Full Size | PDF

 

Fig. 5 Similar to Fig. 4, but for a 6-element cubic array and corresponding homogeneous cube. See also Media 4 for time-domain animations in the y-z plane.

Download Full Size | PDF

 

Fig. 6 Similar to Figs. 4–6, but for a 10-element cubic array and corresponding homogeneous cube. See also Media 5 for time-domain animations in the y-z plane.

Download Full Size | PDF

The equivalent homogeneous cubes were chosen to be of side lengths $3.9d$, $5.9d$, and $9.9d$ for cubes composed of $4^3$, $6^3$, and ${10}^3$ elements, respectively. Obviously, the equivalent homogenized parameters, ${\epsilon }_2$ and ${\mu }_2$, are the same as those used for the metamaterial disk. The same mesh parameters are also used, leading to $N=105984$, $N=357696$, and $N=1656000$ edges for arrays with 4, 6, and 10 elements per side whereas the corresponding homogenized cubes require just $N=2592$, $N=5832$, and $N=16200$ edges, respectively. The computational costs for the various simulations are compared in Table 2, again demonstrating large computational advantages in the homogenized cases.

Table 2. Computational requirements for the simulations in Figs. 4-6.

View Table | View all tables in this article

The calculated RCS patterns and the electric near-field distributions on the two central symmetry planes for each size are shown in panels (a)-(b) and (c)-(h) of Figs. 4–6, respectively. In panels (c)-(h) of Figs. 4–6, the left column refers to the metamaterial array, the center column to the homogenized cube, and the right column shows the normalized field difference. The near-field distributions show good visual agreement in all three planes of polarization. The far-field results also show good agreement in the two main planes. The complexity of the RCS pattern grows with the array size, with multiple side lobes in both planes, and the RCS pattern of the homogenized model predicts the RCS features well, including the minor side lobes.

A visual inspection of the results for the 10³-element cube in Fig. 6 indicates that the homogenization model can predict to a large degree the far-field scattering of large finite-size 3-D metamaterial arrays using bulk equivalent parameters, allowing huge computational savings. Even for smaller cube-shaped arrays, good agreement is observed and the main physics of the problem and the overall near-field effects appear to be qualitatively captured. It appears that the straight edges in this configuration improve the homogenized model performance by avoiding granularity effects and truncations at the edges that were necessary for a curved boundary. It is interesting to note that the scattering patterns for the homogenized cubes show less sharp zeros than the corresponding arrays. Still, the minima and position of the lobes are well aligned in all configurations. Panel (i) in Figs. 4–6 shows a snapshot in the time-domain electric field around and inside the arrays of different size. In the smallest array, a standing wave is observed inside the metamaterial array and homogenized cube due to the mismatch at the two edges. The fields are very well captured right outside the boundaries of the object. As the array becomes larger, negative slope and backward phase propagation become more evident inside the structures, consistent with the double-negative nature of the metamaterial. The near- and far-field distributions are well captured for all sizes of the metamaterial objects, implying that the homogenized model can be used with confidence in modeling large metamaterial arrays.

6. Conclusions

This paper investigated the applicability to finite-size metamaterial arrays of a homogenized model based on dipolar approximations and infinite periodic arrays. The near- and far-field scattering response of large-scale arrays of densely packed magnetodielectric spheres arranged in cylindrical and cubic shapes were numerically calculated using a fast surface integral equation solver at frequencies for which the bulk metamaterial supports positive or negative index of refraction. Comparing scattering patterns and near-field distributions inside and outside the structures, we have shown that the homogenized model can give a reasonably good description of the propagation and scattering properties of finite-size metamaterial objects, even for negative-index, near-resonance operation. Accuracy of the results have been shown to improve for larger arrays, for which edge effects are less important, and for straight edges, for which staircasing effects are not observed at the array boundary. The results show that a rigorous homogenization model can correctly capture the main physics of the scattering from finite metamaterial objects and lead to large computational savings in the modeling, design and application of metamaterials in realistic devices. Although the metamaterial arrays used to verify the homogenization theory in this work are composed of magnetodielectric spheres that are used to realize negative index metamaterials in the microwave regime, our theory can be directly applied to optical and terahertz frequencies, as long as subwavelength inclusions are considered.