Using a discrete dipole approximation to predict complete scattering of complicated metamaterials

Over the past decade, the metamaterial concept has been introduced and has seen rapid adoption by the community. The term 'metamaterial' [1] is a modern one, referring to an artificially structured medium formed by a collection of elements, each of which is significantly smaller than the operating wavelength, and which has properties distinct from the constituents [2, 3]. In the metamaterial limit where the element-to-wavelength ratio is small enough (typically λ/10), the collection of elements can conceptually be replaced by a continuous medium characterized by constitutive parameters.

One of the distinguishing characteristics of recent metamaterials research, relative to earlier effective medium and homogenization theories [4, 5], is the use of full-wave numerical simulations to arrive at constitutive parameters for elements with arbitrary geometry and composed of arbitrary materials. In a commonly used metamaterial design approach, a numerical retrieval is performed to determine the effective constitutive parameters from consideration of the computed reflection and transmission coefficients (or scattering parameters) for a single layer of elements. The computed scattering (S-) parameters for the layer of elements can be related to the Fresnel scattering coefficients for an equivalent thickness of continuous material [6]. The scattering (S-) parameters retrieval method has proven extremely efficient in metamaterial design, since only one unit cell of the planar slab need be meshed in the full-wave solver, with the rest of the planar metamaterial array implicit through the application of periodic boundary conditions. Complicated metamaterial structures can thus be designed by first simulating a hypothetical continuous material using a full-wave solver, and then designing elements by the retrieval approach to approximate the desired structure. In this way, a computationally difficult full-wave simulation including all details of hundreds or thousands of metamaterial elements is avoided.

A key simplifying assumption implicit in the S-parameters retrieval approach (as well as the related field-averaging retrieval approach) is that the interaction between metamaterial elements is predominantly dipolar. While it is not difficult to violate this assumption in an arbitrary metamaterial, the validity of the dipole approximation is critical for the successful application of common retrieval methods. When the higher order multipoles of nearby elements begin to couple significantly, the effective constitutive parameters determined from S-parameter as well as many other retrieval methods will fail to accurately describe the composite metamaterial.

Thus, although the S-parameter retrieval approach can, at best, be considered approximate, nevertheless it has been utilized with great success in the design of complex metamaterial devices, and serves as a useful starting point to harness the enormous palette of material response available from metamaterial elements. While more complicated techniques can provide a rigorous description of scattering behavior [7–9], the application of such techniques tends to obviate the advantage of the metamaterial approach—which is the simplification of the design of complex electromagnetic media. Fortunately, in practice it is generally not difficult to realize elements with predominately dipolar interaction, as has been demonstrated in an increasing number of metamaterial experiments.

The numerical homogenization-based retrieval techniques (such as S-parameter retrieval) discussed above are useful because the details of the dipolar response and interaction of elements within a periodic medium are subsumed into the resulting effective constitutive parameters. However, this approach also carries with it certain inconvenient assumptions implicit within the retrieval method: including a global-symmetry enforced by the periodic boundaries and an assumption that there is no magnetoelectric coupling between metamaterial elements. The validity of such retrieval methods can thus become suspect when gradient or transformation optical media are being designed, where the local environment of each element does not necessarily correspond to any crystalline structure. While retrievals performed with the assumption of periodicity have been successfully employed in the design of both gradient index [10–12] and transformation optical [13–18] media, the question remains as to the extent of error this retrieval approach introduces. Moreover, when magnetoelectric coupling between elements (or intra-element) is significant, the standard retrieval results provide an inherently flawed description of the medium [19, 20].

It is always possible, in principle, to determine the exact electromagnetic response of a collection of metamaterial elements through a full-wave simulation whose domain includes all of the elements arranged in the configuration of interest. Though full-wave simulations on metamaterials have been reported that include large numbers of elements within the simulation domain [21], such simulations still remain computationally expensive and negate the advantage of simple metamaterial design. The difficulty is even greater for nanophotonic metamaterials where the detailed response of all constituent materials, which are typically frequency dispersive, becomes important and the computational requirements to simulate a large scale metamaterial device become impractical.

Fortunately, there is an alternative route to simulating the response of an arbitrary arrangement of meta-atoms. Given that the dipolar response is the key interaction between metamaterial elements [22], a hybrid modeling approach can be considered in which the electric and magnetic dipole responses associated with an arbitrary metamaterial element are first determined, and then the properties of the metamaterial device comprising the collection of metamaterial elements is then computed self-consistently. The concept of extracting the dipolar response for a metamaterial element has been recently discussed [22, 23] with the goal of achieving a metamaterial description that does not rely on the assumptions tied to the past retrieval approaches.

Once electric and magnetic polarizabilities can be ascribed to metamaterial elements, existing computational methods such as the discrete-dipole approximation (DDA) [24–26] can be applied to determine the field distributions, scattering properties and any other electromagnetic detail of the metamaterial composite. With the correct polarizabilities assigned to each element, metamaterials with tens of thousands of elements can be simulated; and, in principle, the simulated structures can include all aspects of the physical structures—including magnetoelectric coupling and spatial dispersion—to the extent that dipole coupling predominates. For gradient and transformation optical media, the DDA approach avoids questions regarding the validity of retrieved effective medium parameters in the face of asymmetry, and provides an alternative approach for the global optimization of complex metamaterial media.

In this paper, we approach the problem of simulating complex structures of meta-atoms by considering each metamaterial element—or meta-atom—as a localized point-particle with a finite number of multipole terms. Then, we express the fields from each meta-atom in terms of a multipole expansion of electric and magnetic polarizabilities, rather than as a homogenized effective material. An advantage of this approach is that when a collection of meta-atoms is placed in a regular array, DDA can be leveraged to quickly calculate the behavior of the complete metamaterial structure. The geometric layout and unit-cell construction of metamaterials makes the DDA a seemingly ideal tool for metamaterials.

The reduction of meta-atom to dipole moment is illustrated in figure 1. In figure 1(a), a meta-atom is shown in the absence of any other interacting elements. An inherent polarizability can be determined for the meta-atom, for example, by assuming open boundary conditions and calculating the scattered dipole strength under plane wave illumination. More common is the situation shown in figure 1(b), in which the meta-atom is infinitely repeated in a cubic lattice. Because the total field acting on a given meta-atom includes the interaction fields from all other responding elements, the extraction of the polarizability must take into account an interaction factor. This step is naturally taken in numerical retrievals, but the global symmetry remains an inherent assumption. The final goal, however, is that each metamaterial element be conceptually replaced by a point dipole scatterer, as indicated in figure 1(c).

Zoom In Zoom Out Reset image size

Many meta-atoms (including the ubiquitous split-ring resonator (SRR)) possess both an electric and a magnetic response, and we thus accordingly include both electric and magnetic dipole terms in our DDA. Though it has long been an assumption in retrieval methods, it has recently been shown explicitly that the dominant term for SRR meta-atoms near resonance is the dipole term [22], and therefore we restrict the number of multipoles in our expansion to the dipole terms only. For the majority of contemporary metamaterial element geometries, the dipole terms should provide good accuracy, but we note that the method could be extended to include higher order multiple terms to yield more accurate results [27].

Because much of the numerical and experimental confirmation of transformation optical and other metamaterial media has been done in planar, two-dimensional (2D) configurations, we find it useful to explore a 2D version of the DDA to maintain consistency with prior studies. The method, developed in section 2, is achieved by reducing the general expressions for the three dimensional Green's function to a 2D case, summing over point dipoles along an infinite vertical column. This approach is distinct from the well-known standard 2D Green's function which would incorrectly represent a continuous and homogeneous vertical column of dipoles. To the best of our knowledge, the development we present represents the first use of the 2D DDA, and offers significant improvements in computational performance over a full 3D solution for 2D structures.

In section 3 we validate our DDA method by computing the scattering cross section (SCS) of a 2D array of identical dipoles arranged to form a cylinder with constant dielectric. We compare this SCS with that computed from Mie theory for a homogeneous dielectric cylinder. We then apply the 2D DDA in the analysis of several metamaterial cloaking structures similar to those which have been experimentally demonstrated [13]. The power and usefulness of the DDA method is demonstrated, as we investigate the importance of the layout of metamaterial elements in the performance of such metamaterial cloaks. To date, the only other numerical tool suited to fully interrogate the effect of the metamaterial layout on scattering is full-wave simulations— a time-consuming approach that can take days (or weeks) per simulation. In contrast, the DDA technique can be performed quickly and multiple metamaterial configurations can be analyzed within minutes.

As discussed above, the core of the metamaterial DDA approach is the representation of each element in the metamaterial as an electric and a magnetic point-dipole. The interaction of these dipolar-elements with an incident field and with each other is expressed as a Green's function summation in 3D. Derivations of the 3D DDA for electromagnetism can be found in many works, for instance Mulholland et al [28]. We replicate one such derivation for completeness, but relegate it to appendix A for conciseness. Starting directly from the 3D summation of the point-dipole Green's functions

$$\begin{equation} (\bar{\alpha}_i^E)^{-1}\vec{p}_i=\vec{E}_0(\vec{r}_i)+\sum\limits_{j\not=i}\skew5\bar{C}(\vec{r}_{ij})\vec{p}_j-\sum\limits_{j\not=i}\skew5\bar{f}(\vec{r}_{ij})\vec{m}_j, \quad \end{equation} \tag{ 1 }$$

$$\begin{equation} (\bar{\alpha}_i^H)^{-1}\vec{m}_i= \vec{H}_0(\vec{r}_i) +\sum\limits_{j\not=i}\skew5\bar{f}(\vec{r}_{ij})\vec{p}_j+\sum\limits_{j\not=i}\skew5\bar{C}(\vec{r}_{ij})\vec{m}_j , \end{equation} \tag{ 2 }$$

the total field experienced by one point dipole is calculated as the incident field plus the fields from all the other dipoles. In most implementations of the DDA, these equations are solved for the electric field at the location of each particle by expressing the polarizations on the left-hand side in terms of the electric and magnetic fields [28, 29]. But in our case, we want to solve for the polarizations themselves and use the set to generate field maps and scattering coefficients.

$$\begin{equation} \left(\left[\begin{array}{@{}ccc@{}} (\alpha^E)^{-1} & 0\\ H 0 & (\alpha^H)^{-1}\\ \end{array}\right] -\left[\begin{array}{@{}ccc@{}} \skew5\bar{C} & -\skew5\bar{f}\\ \skew5\bar{f} & \skew5\bar{C}\\ \end{array}\right]\right) \left[\begin{array}{@{}c@{}} \vec{p}\\ \vec{m}\\ \end{array}\right] =\left[\begin{array}{@{}c@{}} \vec{E}_0\\ \vec{H}_0\\ \end{array}\right]. \end{equation} \tag{ 3 }$$

For our purposes $\bar {C}(0)=\bar {f}(0)=0$, which physically means that the dipole is not affected by its own fields. However, it is worth noting that many realistic meta-atom configurations have a strong local magnetoelectric coupling that could make this term deviate from zero.

In principle, equations (1) and (2) can be solved directly for any 3D collection. However, we wish to consider a subset of metamaterial problems where light propagates in a plane (which we take as the x–y plane), incident upon a collection of dipoles each of which extends infinitely and homogeneously out of the plane (z direction). This seemingly contrived situation is relevant to a wide range of previously reported experiments and simulations, as quasi-2D metamaterials allow for controlled polarization and simple fabrication. In this case, the lattice can be viewed as a finite set of infinite columns of identical dipoles with their axes pointing out of the plane. Experimentally, this is realized by using TE light (electric field vector orthogonal to wavevector) and constraining the sample between two conducting slabs. The electric and magnetic polarizability tensors are left to be completely arbitrary.

For the supposed quasi-2D geometry, our goal is to reduce the matrices $\bar {f}(r)$ and $\bar {C}(r)$ that represent the fields of individual dipoles to the effective fields $\bar {f}^{\mathrm {eff}}(r)$ and $\bar {C}^{\mathrm {eff}}(r)$ which represent an infinite column of dipoles. This is easily done by numerically performing a sum over the column dipoles in the z direction, continuing until a satisfactory level of convergence is reached. The resulting Green's function arrived at by this summation is distinctly different from the standard 2D Green's function found if the sum is replaced by an integral. Having performed this summation, the entire problem can be re-expressed using the effective fields and polarizabilites of infinite dipole columns. We know by symmetry that $\vec {m}$ and $\vec {p}$ must be constant in the z direction, so they can be factored out to the right of the summation in z.

$$\begin{eqnarray} &&\fl (\bar{\alpha}_i^E)^{-1}\vec{p}_i -\sum\limits_{\kappa \not=0}\skew5\bar{C}(\kappa l\hat z)\vec{p}_i -\sum\limits_{j\not=i}\left[\sum\limits_{\kappa}\skew5\bar{C}(\vec{r}_{ij}+\kappa l\hat z)\right]\vec{p}_j +\sum\limits_{j\not=i}\left[\sum\limits_{\kappa }\skew5\bar{f}(\vec{r}_{ij}+\kappa l\hat z)\right]\vec{m}_j\nonumber\\ &&+\sum\limits_{\kappa \not=0}\skew5\bar{f}(\kappa l\hat z)\vec{m}_i =\vec{E}_0(\vec{r}_i), \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &&\fl (\bar{\alpha}_i^H)^{-1}\vec{m}_i -\sum\limits_{\kappa \not=0}\skew5\bar{f}(\kappa l\hat z)\vec{p}_i -\sum\limits_{j\not=i}\left[\sum\limits_\kappa \skew5\bar{f}(\vec{r}_{ij}+\kappa l\hat z)\right]\vec{p}_j -\sum\limits_{j\not=i}\left[\sum\limits_\kappa \skew5\bar{C}(\vec{r}_{ij}+\kappa l\hat z)\right]\vec{m}_j\nonumber\\ &&-\sum\limits_{\kappa \not=0}\skew5\bar{C}(\kappa l\hat z)\vec{m}_i =\vec{H}_0(\vec{r}_i). \end{eqnarray} \tag{ 5 }$$

The indices i and j are now limited to particles in the x–y plane, while the sum over κ accounts for the response above and below the plane, and l is the lattice constant. Note that the last summation in each of equations (4) and (5) is for the column of dipoles directly above and below dipole i. This sum over m describes the interaction constant of a periodic column of dipoles as a function of distance from the column—analogous to the interaction constant between sheets of regularly spaced dipoles [30, 31].

Memory and processing time can be saved by noting at the outset that there is anti-symmetry about the x–y plane in the matrices that will cause some of the elements to go to zero within $\bar f$ and $\bar C$ (see equations (A.15) and (A.16) in appendix A) when they are summed in the z direction.

$$\begin{equation} \skew5\bar{C}(r)=\left[\begin{array}{@{}ccc@{}} a(r)+b(r)\hat r_x^2 & b(r)\hat r_x\hat r_y & 0\\ b(r)\hat r_x\hat r_y & a(r)+b(r)\hat r_y^2 & 0\\ 0 & 0 & a(r)+b(r)\hat r_z^2\end{array}\right], \end{equation} \tag{ 6 }$$

$$\begin{equation} \skew5\bar{f}(r)=d(r)\left[\begin{array}{@{}ccc@{}} 0 & 0 & \hat r_y\\ 0 & 0 & -\hat r_x\\ -\hat r_y & \hat r_x & 0\end{array}\right]. \end{equation} \tag{ 7 }$$

$\bar C$ can be expressed in terms of three matrices that are independent of z, multiplied by functions that are dependent on z.

$$\begin{equation} \bar C(r)=a(r)\bar I+\frac{b(r)}{r^2}\left[\begin{array}{@{}ccc@{}} x^2 & xy & 0\\ xy & y^2 & 0\\ 0 & 0 & 0\end{array}\right] +b(r)\frac{z^2}{r^2}\left[\begin{array}{@{}ccc@{}} 0 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 1\end{array}\right], \end{equation} \tag{ 8 }$$

$$\begin{equation} \skew5\bar{f}(r)=\frac{d(r)}{r}\left[\begin{array}{@{}ccc@{}} 0 & 0 & y\\ 0 & 0 & -x\\ -y & x & 0\end{array}\right]. \end{equation} \tag{ 9 }$$

Now $\bar {f}^{\mathrm {eff}}$ and $\bar {C}^{\mathrm {eff}}$ can be completely constructed by first summing the four z-dependent matrix coefficients in these equations over the z direction until they converge, and then applying equations (A.2) and (A.3). By noting that C_xy = C_yx, f_zx = −f_xz, and f_yz = −f_zy, we can completely reduce $\bar C$ and $\bar f$ to six scalar functions: C_xx,C_xy,C_yy,C_zz,f_xz and f_yz.

The dipole lattice columns $\sum \nolimits _{\kappa \not =0}\bar {C}(\kappa l\hat z)$ can also be accounted for by an effective polarizability for a column of dipoles. Following the symmetry of the system,

$$\begin{equation} \sum\limits_{\kappa \not=0}C_{xx}(\kappa l\hat z)=\sum\limits_{\kappa \not=0}C_{yy}(ml\hat z) =\sum\limits_{\kappa \not=0}a(\kappa l\hat z)=\gamma_x, \end{equation} \tag{ 10 }$$

$$\begin{equation} \sum\limits_{\kappa \not=0}C_{zz}(\kappa l\hat z)=\sum\limits_{\kappa \not=0}C_{zz}(\kappa l\hat z) =\sum\limits_{\kappa \not=0}\big(a(\kappa l\hat z)+b(\kappa l\hat z)\big)=\gamma_z, \end{equation} \tag{ 11 }$$

$$\begin{equation} \sum\limits_{\kappa \not=0}C_{xy}(\kappa l\hat z)=\sum\limits_{\kappa \not=0}C_{yx}(\kappa l\hat z)=0. \end{equation} \tag{ 12 }$$

For a 2D simulation with the electric field in the plane, the sum $\sum \nolimits _{\kappa \not =0}\bar {f}(\kappa l\hat z)$ is identically zero. Once the correction factors are known, the effective polarizability for a column of dipoles can now be expressed as,

$$\begin{equation} (\alpha^{\mathrm{eff}})^{-1}=(\alpha)^{-1}- \left[\begin{array}{@{}ccc@{}} \gamma_x & 0 & 0\\ 0 & \gamma_x & 0\\ 0 & 0 & \gamma_z\\ \end{array}\right]. \end{equation} \tag{ 13 }$$

If we choose to restrict α and β such that the electric polarization is in the x–y plane while the magnetic polarization is in the z direction, then the matrix in equation (13) representing the system is reduced to 3N × 3N. The same kind of expression can be developed using the same method to find a correction for a sheet of dipoles in a square lattice for the 1D problem.

We first validate the 2D DDA approach by simulating a scattering geometry having a known, closed form analytical solution. As a computed metric for comparison, we choose the SCS, which provides a simple and key figure-of-merit applicable also to the more complicated cloaking structures that will be subsequently investigated. The solution of a dielectric cylinder, infinite in extent, can be found from standard Mie theory [32].

In order to apply the DDA method of section 2, it is necessary to assign the correct electric and magnetic polarizabilities (equations (1) and (2)) to the dipole elements that will result in the collective behaving as a uniform dielectric medium. The step of assigning polarizabilities is a subtlety of the DDA method and must be performed cautiously, since the correct choice of polarizability must reflect the influence of all interacting dipoles. An initial approach makes use of the Clausius–Mossotti equation, which has the form

$$\begin{equation} \bar{\alpha}^{\mathrm{CM}}=V(\bar{\epsilon}-I)(\bar{\epsilon}-C+C\bar{\epsilon})^{-1}, \end{equation} \tag{ 14 }$$

where V is the volume around the element, I is the identity matrix, and C is the depolarization or interaction factor, which generally is a tensor. For a dipole located deep within the interior of a cubic lattice of identical dipoles, $C=\frac {1}{3}I$. For a dipole located on one of the outermost surface layers, there are corrections for the interaction factor that take into account the asymmetric environment [30]. Likewise, for lattices having symmetries other than cubic, or for non-periodic dipole arrays, the interaction factors will be different and may generally vary from dipole to dipole.

For a cubic lattice, the Clausius–Mossotti equation is exact in the static limit. However, outside of the static limit, it can be shown that the Clausius–Mossotti equation applied to the DDA results in a collective structure that violates the optical theorem. In most implementations of the DDA, this issue is addressed by setting the imaginary part of the polarizability such that it satisfies conservation of energy [29, 30].

$$\begin{equation} \alpha =\alpha^{\mathrm{CM}}\left(I-\mathrm{i}\frac{k^3}{6\pi}\alpha^{\mathrm{CM}}\right)^{-1}. \end{equation} \tag{ 15 }$$

This correction gives an exact solution for the imaginary part of the polarizability of lossless dipoles in the nonzero frequency regime, but the real part of the polarizability is still an approximation. Various methods exist for improving this approximation, [29, 33] but for our purposes equation (15) is sufficiently accurate.

Now that we can confidently express the homogeneous material parameters of the cylinder in terms of the complex polarizability of a cubic lattice of dipoles, we proceed with calculation of (SCS) for a dielectric cylinder using the DDA. In most traditional DDA implementations, the SCS for arbitrary objects is calculated by first finding the difference of the absorption and extinction of each individual dipole, and then summing the contribution from all the dipoles. However, given that we are simulating an infinite column of dipoles, this method becomes cumbersome. Instead, we calculate the cross section of the dielectric cylinder using a Bessel function expansion of the scattered fields computed from the DDA approximation to the cylinder. For an incident plane wave with its magnetic field oriented along the cylinder axis, the scattering becomes a scalar problem in terms of the magnetic field [32].

$$\begin{equation} H_s(r,\phi)=|H_0|\sum\limits_{\kappa =-\infty}^\infty a_{\kappa}H_{\kappa}^1(kr)\,\mathrm{e}^{-\mathrm{i}\kappa\phi}, \end{equation} \tag{ 16 }$$

$$\begin{equation} a_{\kappa}=\frac{1}{2\pi |H_0|H_{\kappa}^1(kr)}\int\nolimits_0^{2\pi}H_s\,\mathrm{e}^{-\mathrm{i}\kappa\phi}, \end{equation} \tag{ 17 }$$

$$\begin{equation} C_{\mathrm{sca}}=\frac{4}{k|H_0|^2}\sum\limits_{\kappa=-\infty}^\infty |a_{\kappa}|^2. \end{equation} \tag{ 18 }$$

Figure 2(a) depicts the scattering cylinder comprising of a collection of dipoles. We select a cubic arrangement of dipoles, with lattice constant 0.09λ, and with moments weighted to have an effective refractive index of n = 1.5. Figure 2(b) shows the field pattern derived from the DDA for a cylinder of radius 1.35λ. Figure 2(c) (red) shows the calculated SCS using equation (18) for various cylinder radii alongside analytical results from Mie theory. In this geometry, SCS has the units of area (meters). We also re-run the scattering at a denser lattice of 0.045λ (green), and see subtle improvement but essentially the same SCS. The results indicate that 0.09λ is sufficiently sub-wavelength to form a good effective medium (at these polarizabilities) for the scattering cylinder. Overall the observed good agreement for the SCS between the 2D DDA and Mie theory validates our approach.

Zoom In Zoom Out Reset image size

Having validated the accuracy of our DDA against a known test case, we now wish to tackle more complex media. The electromagnetic cloak [34] is one of the more captivating designs to come out of the field of transformations optics, and implementations of electromagnetic cloaks are among the most complicated metamaterial devices [13] realized to date. The total SCS is an appropriate metric for the performance of a cloak, and has been previously used to benchmark both measured and simulated cloak variations [35, 36]. To illustrate the 2D DDA method, we replicate the parameters for the cloak design experimentally demonstrated in [13]. Initially we consider cloaks with vacuum in the interior, such that the SCS in the absence of the cloak would be zero.

In a manner similar to that done for the scattering cylinder in figure 2, we use the Clausius–Mossotti equation to determine the polarizabilities for each of the meta-atoms used to form the cloak [13]. Unlike the cylinder, however, which was assumed to be a uniform dielectric material, the cloak is a transformation optical medium for which the constitutive tensor elements generally vary spatially with radius and according to the form of the cloak.

$$\begin{equation} \epsilon_r = \mu_r = \frac{r-a}{r} , \end{equation} \tag{ 19 }$$

$$\begin{equation} \epsilon_\theta = \mu_\theta = \frac{r}{r-a}, \end{equation} \tag{ 20 }$$

$$\begin{equation}\epsilon_z = \mu_z = \left(\frac{b}{b-a}\right) ^2 \frac{r-a}{r}.\end{equation} \tag{ 21 }$$

We must therefore determine the polarizabilities for each dipole at every position throughout the cloak in order to apply the DDA, in the absence of any periodicity. That is, each dipole is surrounded by dipoles with distinct polarizabilities. Since the polarizabilities and Green's functions in equations (4) and (5) are expressed in rectangular coordinates, the material parameters for the cloak must be transformed into the rectangular basis before calculating the polarizabilities from the Clausius–Mossotti equation.

The singularity in equation (20) was handled for each of our cloaks by removing any dipoles that lie along the circle r = a. This is an approximation that must be made for any physical metamaterial cloak that uses this transformation design, and so the approximation is one that would realistically predict the necessary scattering properties of a cloak made with finite materials.

To easily use the DDA, the polarizability for each dipole must first be calculated using the Clausius–Mossotti equation with an interaction factor determined by implementing certain assumptions. As a first set of assumptions, we assume the dipoles are laid out in a cubic lattice, and utilize the interaction factor assuming a given dipole is arranged in an infinite cubic lattice. This assumption is equivalent to that used in S-parameter retrievals (such as was done for the cloak in [13]), where the assumption of perfect infinite cubic-lattices is enforced via periodic boundary conditions used in full wave solvers. By using the Clausius–Mossotti equation with an interaction factor of 1/3, we ensure that homogenization assumptions on the unit-cell or meta-atom level are the same between the S-parameter and DDA methods. If we then rearrange the meta-atoms in a simulation, but use the same interaction factors, we introduce the same error that currently occurs in metamaterial structures designed using S-parameter retrieval. Thus, the DDA allows us to investigate the errors caused by various methods of positioning meta-atoms in a unique manner. Figure 3(a) depicts a cubic-lattice metamaterial arranged into a 2D cloak. We also plot the scattered fields (planewave incident) and tabulate the total SCS_cubic = 0.598.

Zoom In Zoom Out Reset image size

Due to the difficulty in physically orienting unit-cells axisymmetrically when placed on a Cartesian grid (such as in figure 3(a), this cubic design has not been utilized for experimental cloaks [13]. Instead, meta-atoms have been arranged along circumferential rings. Since the circumferential arrangement deviates from the S-parameter-retrieval geometry, it was previously unknown how the modified inter-dipole coupling would affect the scattering performance. The hope during S-parameter design has usually been that if the arrangement is nearly cubic on a local scale, the fabricated non-cubic structures may still follow the cubic retrieval values. Using the DDA method we can investigate the impact of locally and globally non-cubic lattices. Figure 3(b) depicts a circumferential dipole arrangement like that fabricated in [13], and plots the scattered fields. Calculating the scattering, we find SCS_circ = 0.705, slightly higher than that of the cubic layout. Although this circumferential arrangement lacks the exact cubic-layout of (a), it has an axisymmetric symmetry that likely maintains some measure of crystallinity locally. From a comparison of the two SCS values, it is apparent that the layout of the meta-atoms has a significant effect on the scattering of the composite, presumably because the interaction factors for the dipoles have deviated significantly from the assumed value of 1/3. One other meta-atom arrangement is of particular interest, due to its ease of fabrication: the case where meta-atoms are located along radial arms. In this arrangement, meta-atoms are still spaced evenly outward in radius. However, as the total number of meta-atoms remains constant, the inter-atom spacing varies circumferentially as $\frac {1}{r}$. This case is depicted in figure 3(c), where the scattered fields are also shown. This case turns out to be the worst of the three lattices, with SCS_rad = 0.924. Since this arrangement is the least crystalline, its relatively poor performance is not surprising. However, as we will see below, the cloak formed from dipoles laid out in a radial pattern still reduces the net SCS of an object contained within the cloaking volume; the relative fabrication ease of this design may still enable it to be useful in some situations.

Finally, we consider the three lattice geometries when a scattering center is included within the cloaking volume. For consistency we use the cubic-layout dielectric (n = 1.5) cylinder from figure 2, choosing a radius of 1.35. The performance metric of the cloak now becomes the reduction of the cross-section below SCS_scatt = 5.35. Figure 4 shows the performance of the three lattice types for this cloak. As before, the cubic-lattice cloak performs the best (SCS_cubic+scatt = 0.669) followed by the circumferential cloak (SCS_circ+scatt = 0.830), with the radial cloak performing worst (SCS_rad+scatt = 1.175). The fact that the cubic-layout shows a SCS 20% lower than the geometry used in experiments [13] suggests there is significant room for optimization of metamaterial devices by adjusting the layout of unit-cells.

Zoom In Zoom Out Reset image size

As always, the unique meso-scale construction of artificially structured metamaterials adds design complications that would be absent were conventional materials available. Spatial dispersion, for example, is almost always present in metamaterials, and is a major consideration in their simulation and design. The S-parameters retrieval is known to be incorrect when applied to gradient and transformation optical media, as it incorrectly assumes global symmetries and artificially enforces homogenization. The DDA tool presented in this paper allows us to escape some of these problems when simulating the operation of a full metamaterial device—especially those problems related to lattice symmetry and inter meta-atom coupling. These effects are important considerations for metamaterial device design; change of lattice-geometry within the same effective medium prescription has been shown here to correspond to almost a 100% change in the SCS of an electromagnetic cloak.

A key point in applying the DDA to simulate a metamaterial structure is that the polarizabilities for each of the dipolar elements must first be determined. Though we have shown the Clausius–Mossotti equation is a good starting point, and validates the S-parameters-based design approach that has been applied, it is still clear that errors exist with this method. In principle, there is no reason for the performance of the three cloaks simulated above to differ, were the correct polarizabilities found for each of the dipoles in the particular layouts. As a future refinement of the DDA approach, it should be possible to start a simulation with the polarizabilities of each dipole determined from the Clausius–Mossotti equation, and then optimize the interaction factors while minimizing the SCS or other desired figure-of-merit. This global optimization approach should lead to transformation optical media with much greater performance characteristics.

Despite the advantages of the DDA, complications nevertheless exist. At the heart of these problems is that in our employment of Green's function we have assumed magnetoelectric point-dipoles. For many meta-atom designs, the existence of geometric asymmetry and self-coupling makes this incorrect. The response of the metamaterial element is spatially distributed, and in the presence of certain asymmetries cannot be condensed to a response at a single point, even if all multipole moments are allowed. Along these lines, higher order modes (beyond dipole) in the multipole expansion can be significant in some meta-atoms, and these are not taken into account by DDA. Methods exist to address such concerns. For instance, interdigitation of several offset lattices of dipoles can reproduce the self-coupling of meta-atom geometries. Similarly, there exist DDA-like methods that include higher order terms in the multipole expansion [27]. However, when considering such improvements we must keep in mind that the full solution is always available via (painfully slow) full-wave solvers—and what we are after is a balance between speed and accuracy. Realistic design of metamaterial devices for application requires a fast feedback loop between effective medium prescription, meta-atom layout, and meta-atom geometry, allowing an iterative process to optimize towards a performance metric. Full-wave solvers have traditionally proved too cumbersome to include in such a feedback loop, and while the S-parameter retrieval has served surprisingly well to date, it cannot address meta-atom layout considerations. It is the second of these tri-scale considerations necessary for metamaterial design (effective medium, lattice, meta-atom) that the DDA tool is uniquely situated to address: the meta-atom layout, and characterizing the way in which interactions between meta-atoms affect the performance of the composite.

This work was supported through a Multidisciplinary University Research Initiative, sponsored by the Army Research Office (grant no. W911NF-09-1-0539). TD acknowledges support from an IC postdoctoral fellowship.

In this appendix, the expression for the 3D DDA is rederived in Gaussian units. Following the points discussed in the introduction, we assume point-particles and write the electric and magnetic dipoles ($\vec {p}_i$ and $\vec {m}_i$) in terms of the familiar electric and magnetic polarizability tensors $\bar {\alpha }_E$ and $\bar {\alpha }_H$.

$$\begin{equation} \vec{p}_i=\bar{\alpha}_i^E\vec{E}_i, \end{equation} \tag{ A.1 }$$

$$\begin{equation} \vec{m}_i=\bar{\alpha}_i^H\vec{H}_i. \end{equation}\noindent \tag{ A.2 }$$

The fields for an electric dipole can be expressed as

$$\begin{equation} \vec{E}_p(\vec{r})=\frac{\mathrm{e}^{\mathrm{i}kr}}{r}\Big[\vec{p}\left(k^2-\frac{1}{r^2}+\frac{\mathrm{i}k}{r}\right)+\hat{r}(\hat{r} \cdot \vec{p})\left(-k^2+\frac{3}{r^2}-\frac{\mathrm{i}3k}{r}\right)\Big],\end{equation}\noindent \tag{ A.3 }$$

$$\begin{equation} \vec{H}_p(\vec{r})=(\hat{r}\times\vec{p})\frac{\mathrm{e}^{\mathrm{i}kr}}{r}\left(k^2+\frac{\mathrm{i}k}{r}\right), \end{equation}\noindent \tag{ A.4 }$$

where $\vec {r}$ is a radius vector in 3D space, $\hat {r}$ is it's unit-vector, and $\vec {p}$ is the dipole polarizaton vector. Similarly, the fields for a magnetic dipole can be expressed as

$$\begin{equation} \vec{H}_m(\vec{r})=\frac{\mathrm{e}^{\mathrm{i}kr}}{r}\Bigg[\vec{m}\left(k^2-\frac{1}{r^2}+\frac{\mathrm{i}k}{r}\right)+\hat{r}(\hat{r} \cdot \vec{m})\left(-k^2+\frac{3}{r^2}-\frac{\mathrm{i}3k}{r}\right)\Bigg], \end{equation}\noindent \tag{ A.5 }$$

$$\begin{equation}\vec{E}_m(\vec{r})=-(\hat{r}\times\vec{m})\frac{\mathrm{e}^{\mathrm{i}kr}}{r}\left(k^2+\frac{\mathrm{i}k}{r}\right). \end{equation}\noindent \tag{ A.6 }$$

Following Mulholland et al [28], we can define the functions a(r),b(r),d(r) that depend only on the magnitude of r.

$$\begin{equation} a(r)=\frac{\mathrm{e}^{\mathrm{i}kr}}{r}\left(k^2-\frac{1}{r^2}+\frac{\mathrm{i}k}{r}\right), \end{equation}\noindent \tag{ A.7 }$$

$$\begin{equation} b(r)=\frac{\mathrm{e}^{\mathrm{i}kr}}{r}\left(-k^2+\frac{3}{r^2}-\frac{\mathrm{i}3k}{r}\right),\end{equation}\noindent \tag{ A.8 }$$

$$\begin{equation} d(r)=\frac{\mathrm{e}^{\mathrm{i}kr}}{r}\left(k^2+\frac{\mathrm{i}k}{r}\right). \end{equation} \tag{ A.9 }$$

Now the fields experienced by particle i due to particle j can be expressed as:

$$\begin{equation} (\bar{\alpha}_i^E)^{-1}\vec{p}_i=a(r_{ij})\vec{p_j}+b(r_{ij})\hat{r}_{ij}(\hat{r}_{ij}\cdot\vec{p}_j)-d(r_{ij})(\hat{r}_{ij}\times\vec{m}_j), \end{equation} \noindent \tag{ A.10 }$$

$$\begin{equation} (\bar{\alpha}_i^H)^{-1}\vec{m}_i=a(r_{ij})\vec{m_j}+b(r_{ij})\hat{r}_{ij}(\hat{r}_{ij}\cdot\vec{m}_j)+d(r_{ij})(\hat{r}_{ij}\times\vec{p}_j), \end{equation} \tag{ A.11 }$$

where

$$\begin{equation} r_{ij}=\vec{r}_i-\vec{r}_j. \end{equation} \tag{ A.12 }$$

The vectorial operations and radial functions can all be combined into two matrices that represent the Green's function for this system.

$$\begin{equation} \skew5\bar{C}(r)\vec{p}=a(r)\vec{p}+b(r)\hat{r}(\hat{r}\cdot\vec{p}), \end{equation} \tag{ A.13 }$$

$$\begin{equation} \skew5\bar{f}(r)\vec{p}=d(r)(\hat{r}\times\vec{p}), \end{equation} \tag{ A.14 }$$

$$\begin{equation} \skew5\bar{C}(r)=\left[\begin{array}{@{}ccc@{}} a(r)+b(r)\hat r_x^2 & b(r)\hat r_x\hat r_y & b(r)\hat r_x\hat r_z\\ b(r)\hat r_x\hat r_y & a(r)+b(r)\hat r_y^2 & b(r)\hat r_y\hat r_z\\ b(r)\hat r_z\hat r_x & b(r)\hat r_z\hat r_y & a(r)+b(r)\hat r_z^2\end{array}\right] , \end{equation} \tag{ A.15 }$$

$$\begin{equation} \skew5\bar{f}(r)=d(r)\left[\begin{array}{@{}ccc@{}} 0 & -\hat r_z & \hat r_y\\ \hat r_z & 0 & -\hat r_x\\ -\hat r_y & \hat r_x & 0\end{array}\right] . \end{equation} \noindent \tag{ A.16 }$$

Using these matrix representations for the fields of a point dipole, the total field experienced by one point dipole is calculated as the incident field plus the fields from all the other dipoles.

$$\begin{equation} (\bar{\alpha}_i^E)^{-1}\vec{p}_i=\vec{E}_0(\vec{r}_i)+\sum\limits_{j\not=i}\skew5\bar{C}(\vec{r}_{ij})\vec{p}_j-\sum\limits_{j\not=i}\skew5\bar{f}(\vec{r}_{ij})\vec{m}_j, \end{equation} \tag{ A.17 }$$

$$\begin{equation} (\bar{\alpha}_i^H)^{-1}\vec{m}_i= \vec{H}_0(\vec{r}_i) +\sum\limits_{j\not=i}\skew5\bar{f}(\vec{r}_{ij})\vec{p}_j+\sum\limits_{j\not=i}\skew5\bar{C}(\vec{r}_{ij})\vec{m}_j. \end{equation} \noindent \tag{ A.18 }$$

These two equations comprise the DDA for a 3D collection of scatterers.

The convergence of a discretely summed Green's function can be studied directly, as the function approaches a steady state value. However, the convergence of Green's function is known to be notoriously difficult. Tretyakov's analytical treatment of the problem of summing it in two dimensions involved introducing loss into the system before the integral would converge, and then taking the limit as the loss goes to zero [30]. We show here that the functions do converge without the introduction of loss for a one-dimensional column of dipoles by plotting the real part of the four z-dependent matrix coefficients in equations (8) and (9) for different values of the radial coordinate $\rho =\sqrt {x^2+y^2}$. The convergence of two exemplary functions are given here, and we omit the other two for brevity, noting that the other two converge in a similar fashion.

Zoom In Zoom Out Reset image size