A surface integral equation (SIE) formulation is applied to the analysis of electromagnetic problems involving three-dimensional (3D) piecewise homogenized left-handed metamaterials (LHM). The resulting integral equations are discretized by the well-known method of moments (MoM) and solved via an iterative process. The unknowns are defined only on the interfaces between different media, avoiding the discretization of volumes and surrounding space, which entails a drastic reduction in the number of unknowns arising in the numerical discretization of the equations. Besides, the SIE-MoM formulation inherently includes the radiation condition at infinity, so it is not necessary to artificially include termination absorbing boundary conditions. Some 3D numerical examples are presented to confirm the validity and versatility of this approach on dealing with LHM, also providing some intuitive verifications of the singular properties of these amazing materials.
©2010 Optical Society of America
Double negative materials, with simultaneously negative electric permittivity (ε) and magnetic permeability (μ) form a group of man-made materials that possess non-naturally occurring electromagnetic behavior. They exhibit negative refractive index, supporting backward traveling-waves where the phase propagation is opposite to the flow of energy. In such materials, the electric, magnetic, and propagation vectors form a left-handed system; for this reason they are usually named left-handed materials (LHM), in opposition to the natural right-handed materials (RHM). Some other interesting electrodynamics properties of LHM are the reversal of Snell-Descartes law (by bending the refracted fields in the opposite direction to normal materials), or the reversal of the Doppler effect. Although LHM were theoretically proposed as early as 1968 , they have not been demonstrated until 2000 . In that year, an artificially structured material (or metamaterial) was synthesized and the negative refraction was experimentally demonstrated at microwave frequencies, yielding to the surprising results predicted 30 years before: the material bended electromagnetic fields in the opposite (wrong) direction. Since then, considerable attention has been arisen and remarkable progress has been achieved in this field, leading to the development of potential applications, including mobile communications, medical imaging, “perfect” (superesolution) lenses, cloaking, designing of antennas, and many others. The reader is referred to  for a further discussion and more references of LHM applications.
A typical metamaterial exhibiting negative index is an artificial structure composed of an arrangement of scattering elements (for example metallic wires and split resonators) that are small compared to the wavelength of interest. Such a composite structure provides simultaneously negative permittivity and permeability in a frequency or in a range of frequencies. Due to its structured composition, realistic numerical simulations of metamaterials can be performed by analyzing the actual lattice of microscopic unit elements. In this context, some detailed three-dimensional (3D) modeling of metamaterials has been done using both the differential [4–7] and integral-equation based formulations [8–11]. Nonetheless, in order to observe macroscopic effects, these simulations require the analysis of electrically large objects (with several wavelengths in size), thus composed of a large number of sub-wavelength elements. In fact, to be considered as a medium with collective behavior, the number of unit elements in a 3D arrangement of metamaterial must be the order of thousands per cubed wavelength . All those elements must be modeled accurately in order to obtain the electromagnetic response of the complete structure, which implies the resolution of very large matrix systems. Besides, in the case of integral-equation formulations one has to face the ill-conditioning and slow convergence of the iterative solutions, due to the resonating nature of the microscopic elements . For all the above reasons, even considering the most recent advances both in algorithmic and computer technology, detailed numerical simulations are still restricted to small pieces of metamaterial.
On the other hand, even though metamaterial structures are intrinsically inhomogeneous media, they exhibit bulk properties and can be characterized by effective macroscopic constitutive parameters. There have been reported several homogenization techniques to extract the macroscopic parameters that have proven to be feasible and robust (see for example [13–15]). Consequently, it is very useful, both from the analysis and the design points of view, to study homogenized (continuous) LHM in which the constitutive parameters have been retrieved using homogenization algorithms. This level of abstraction is much more suitable to numerical analysis than the detailed multi-scale approaches, drastically decreasing the amount of required computational resources, and thus allowing the simulation of large structures of metamaterials.
The most common way to analyze homogeneous metamaterials has been the use of differential-equation formulations. Using the finite difference time domain FDTD method, the LHM is usually described in terms of dispersive Drude or Lorentz models specifying the frequency dependence of the material parameters [3,16–20]. In other cases, the finite element method (FEM) was applied to obtain the scattering calculation at a single frequency, avoiding the specification of frequency dependences [21,22]. In any case, the main advantage of the differential-equation techniques is that they can be straightforwardly implemented. However, it is well known that they are volumetric formulations strongly burdened with the discretization of the object and the surrounding space. Besides, it is rather complicated the formulation of efficient absorbers and resistive terminations for negative index media. The above factors limit the applicability of differential methods in the case of large and open radiation and scattering problems.
Nonetheless, a variety of alternative approaches are available for rigorously modeling the interaction of electromagnetic radiation with matter. Concretely, in the case of homogeneous or piecewise homogeneous dielectric objects (such as the homogenized LHM concerned in this work) surface integral-equation (SIE) formulations would be preferable to the differential techniques. With the SIE formulations the problem can be formulated only over the LHM surfaces and interfaces, avoiding the discretization of volumes, and thus providing a drastic reduction in the number of unknowns. Besides, SIE formulations inherently include the radiation condition at infinity, so there is no need at all for termination absorbing boundary conditions. This enables the straightforward analysis of arbitrary radiation and scattering problems in open spaces, which is of great interest in many practical applications of LHM.
In this paper, we apply a SIE formulation, namely the well known method of moments (MoM) , to the simulation of electromagnetic problems involving homogeneous LHM. This approach bears important advantages, especially when the analysis of three-dimensional metamaterial objects in unbounded media becomes of interest. It also brings the possibility of applying the latest breakthrough developments in fast integral-equation methods, such as those presented in [24–26], to the solution of large-scale problems involving LHM. To the best of our knowledge, this is the first time that the integral-equation formulation is adapted to the solution of homogeneous LHM. Additionally, we present simulations in three-dimensions to demonstrate the validity and applicability of the proposed approach. We have focused on problems that provide intuitive verifications of LHM: the Snell-Descartes law experiment, the interaction of a divergent Gaussian beam with a slab of negative index metamaterial, and the interaction of a flat Gaussian beam with a LHM lens. It is worth emphasizing that all the results have been naturally obtained, without enforce any a priori conditions for the propagation inside the metamaterial, which indeed constitutes an additional proof of the singular electromagnetic properties achieved inside LHM.
2. Surface integral-equation formulations of left-handed materials
Surface integral equation methods are usually applied to the solution of three-dimensional electromagnetic scattering problems involving homogeneous or piecewise homogeneous dielectrics of arbitrary shape [27–32]. Integral equations are obtained using the equivalence theorems for closed surfaces delimiting the homogeneous regions. With these techniques, the surfaces and interfaces of the dielectric objects are discretized and the equivalent surface electric and magnetic currents are found in terms of a sum of known vector basis functions. The procedure is described next for a single homogeneous dielectric object. This simplifies the notation without a lack of generality; the procedure is valid for any combination of piecewise composite dielectric objects.
Let us consider a homogeneous dielectric material, and let us denote with the exterior region of the material and with the interior. Let us define with and the constitutive parameters of region, with l = 1, 2. The surface of the object is defined as S, and is the unit vector normal to S pointing toward the interior of region. The electric field integral equation (EFIE) and the magnetic field integral equation (MFIE) can be formulated in two alternative ways, depending on the method applied to project the fields onto the surface S. This leads to the tangential (T) and normal (N) formulations. The tangential integral equations in region are given by:Eq. (6), is the wavenumber in , and we define:
Taking into account the continuity of tangential field components on the surface S, the surface current densities in both sides of S satisfy:
Different known formulations can be obtained depending on the selection of the complex combination parameters , , and [29–32]. For this work, we have chosen , which leads to the JMCFIE combined formulation . JMCFIE is free of internal resonances and it provides a higher level of accuracy for the analysis of LHM objects.
At this point, using the standard MoM procedure the unknown current densities J and M in Eq. (13) are expanded into a sum of known basis functions. The unknown coefficients of these expansions are then obtained by solving a dense matrix system of linear equations, which is derived from the above set of integral equations by applying the Galerkin’s testing procedure (consisting of using the same testing functions as basis functions). Once the solution for the equivalent currents is found, any other parameters of interest (e.g. the near and far field distributions) can be straightforwardly obtained. For the discretization of the problem, in this work we have employed the extensively used Rao-Wilton-Glisson (RWG) basis functions . On the other hand, accurate calculation of the integrals involved in operators of Eq. (5) and Eq. (6) for the RWG functions is crucial for the precision of the method. We have applied Gaussian quadrature rules for the numerical integration of smooth varying integrands, together with the analytical extraction procedures of [34–37] for the accurate evaluation of singular integrals.
Finally, the application of the above-described formulation to the analysis of LHM requires taking special care to the definition of the electromagnetic parameters and properties of the material. Concretely, from the detailed mathematical studies carried out in , if region defines a left-handed medium the following definitions apply for the wavenumber and the intrinsic impedance:
It is worth noting that the previous definitions are valid both for LHM and RHM media. However, the usual definitions for RHM are not valid in the case of LHM.
3. Simulation results
In order to show the validity and applicability of the proposed formulation, we have considered in this section three examples usually found in the literature to provide intuitive verifications of the extraordinary properties of LHM. All of them have been solved on a computer with two 2.5 GHz Intel Xeon 5420 Quad Core processors and 16 GB of memory. The code has been parallelized with OpenMP to take benefit from this multi-core architecture.
The first example is the simulation of the Snell-Descartes law experiment, whose aim is to demonstrate the negative angle of refraction. The usual approach is to determine this refracted angle for an incident wave interacting with a prism comprised of the LHM under test. Usually, a 2D geometry is considered to reduce the computational requirements using FDTD or FEM. In this work, an unrestricted 3D simulation of this experiment has been carried out using the proposed SIE-MoM formulation. A 3D prism has been considered being 10λ0 high, 10λ0 wide, and 4λ0 deep in the greatest dimension, with λ0 the wavelength of the surrounding free-space medium. The angle of the second surface is 8° (this choice was arbitrary). An almost flat Gaussian beam, with a total angular spread of 25° and the waist at a distance of 5λ0 away is orthogonally impinging onto the first surface of the prism from the right hand side. The incident electric field is horizontally polarized.
Figure 1(a) shows the computed electric field intensity for the experiment described above considering a refraction index of n = −1 for the LHM comprised prism. The problem was solved with the SIE-MoM formulation using 40,770 unknowns for the equivalent electric and magnetic currents induced on the surfaces of the prism. The solution time was less than 7 minutes. Once the currents have been calculated, all the other parameters of interest, such as near or far field distributions, can be straightforwardly obtained from them. Looking at the field distribution of Fig. 1(a), the negative angle of refraction due to the LHM prism is clearly observed: the beam is bended in the opposite direction (that is, in the same side as the incident wave with respect to the surface normal). Besides, no reflection effects are observed, since the selected LHM is impedance matched to free-space. Figure 2 shows the angular field pattern of the refracted wave at a distance of 20λ0 away from the prism (in red solid line for the LHM case). It can be observed that the peak of the refracted field occurs at an angle of 196°, which corresponds with 8° from the surface normal. This is exactly equal (and opposite) to the angle of incidence on this surface, as expected from the Snell-Descartes law for n = −1 (the prediction of this law is shown in Fig. 2 as a dotted red straight line). On the other hand, the time evolution of the field inside the LHM prism (not shown in this example) demonstrates that the phase propagation is reversed and that it is in the opposite direction to the power flow (the effect of the phase velocity reversal is shown in the next example of Fig. 3 ).
For the sake of comparison, Fig. 1(b) shows the same above experiment repeated for a prism made up of conventional RHM (we have considered Teflon, with εr = 2.2 and μr = 1). In this case, the refraction is obtained in the positive (conventional) direction. The peak of the refracted beam in Fig. 2 (blue dashed line) occurs at an angle of 176.834°, corresponding to 11.166° from the surface normal. This is in close agreement with the 11.913° predicted by the Snell-Descartes law for the positive index of n = 1.483 (this prediction is shown in Fig. 2 as a dotted blue straight line). Also, note that in this case the prism is not impedance matched to the external medium, hence leading to the usual artifacts due to the wave reflections on the surfaces.
As a second example we have considered a strongly diverging Gaussian beam (with a total angular spread far from the waist of 70°) impinging onto a planar LHM matched 3D slab with n = −1. The slab is 10λ0 high, 10λ0 wide and 2λ0 deep, and the Gaussian beam is normally incident on it from the right hand side. The example was solved in about 8 minutes, using 42,000 unknowns for the electric and magnetic currents induced on the slab surfaces. Figure 3 shows the computed electric field intensity of the beam interacting with the slab. It can be clearly seen how the slab bends the diverging wave towards the beam axis, thus behaving like a lens and focusing the 3D beam (the idea of using planar slabs of LHM to construct lens based on the negative refraction was first suggested by Pendry ). On the other hand, the time-movie of the electric field clearly demonstrates that the phase propagation inside the LHM slab is reversed. It is opposite to the power flow direction. This is a unique and surprising characteristic of LHM.
A planar slab was able to focus a divergent beam in the previous example. However, in order to focus a flat beam (one with nearly an infinite radius of curvature) a curved lens must be applied. Using LHM lens, in contrast to focusing flat beams with a convex lens composed of an RHM medium, we must consider focusing flat beams with a concave lens, because of the reversal of the Snell-Descartes law. So, in the last example we have concerned with a 3D plane-concave lens with n = −1. The lens has been formed by removing a paraboloid volume from the backside of a slab that was 3λ0 deep, 10λ0 high and 10λ0 wide, λ0 being the wavelength of the surrounding free-space medium. The full width of the removed parabolic section at the back face was 6.5λ0. An incident flat Gaussian beam with 15° of total angular spread and the waist at a distance 5λ0 away from the lens is impinging onto its planar side.
The problem was solved using 46,314 unknowns for the electric and magnetic surface currents. The solution time was about 12 minutes. Figure 4 shows the predicted electric field intensity distribution in the horizontal and vertical semi-planes when the intensity is maximum at the focal point (the other semi-planes are symmetrical and they are omitted from the representation for the sake of clarity). It can be seen how the LHM plane-concave lens focuses the beam, contrary to what would happen if it was made of a conventional RHM. On the other hand, it is worth noting that even though the focal point is in the extreme near-field of the lens (due to the lens design and the index of refraction) the focal region is nearly symmetrical in both the vertical and the horizontal plane.
Finally, it is worth mentioning again that the previous results have been obtained directly from Maxwell’s equations and boundary conditions, without assuming any a priory conditions in the formulation of the problem or the propagation inside the materials. Consequently, they constitute a further verification of the fundamental properties of LHM.
A surface integral equation formulation obtained from the Maxwell’s equations and the boundary conditions has been applied to the electromagnetic simulation of problems involving 3D objects comprised of continuous homogenized left-handed metamaterials. The problem was formulated using the JMCFIE formulation. The obtained integral equations were discretized with the MoM, where the unknown equivalent electric and magnetic currents induced on the LHM surfaces and interfaces were represented using RWG vector basis functions. The Galerkin’s method was applied for the testing procedure. Since the SIE-MoM approach is applied only on the surfaces of the objects, it avoids the discretization of volumes that burdens the solution throughout the more usual approaches based on differential-equation formulations. Besides, SIE-MoM formulation inherently includes the radiation condition at infinity, without need of termination conditions. This enables the straightforward analysis of radiation and scattering problems of LHM in unbounded media.
Numerical examples were presented to confirm the validity and versatility of the integral-equation formulation for the resolution of LHM problems. Concretely, three examples usually found in the literature to provide intuitive verifications of the extraordinary properties of these man-made materials were considered: the Snell-Descartes law experiment to demonstrate the negative refraction, the focusing effect of a slab comprised of LHM, and the focusing effect of a plane-concave LHM lens. The simulations were unrestricted in the three dimensions, and they were carried out using much fewer unknowns than it would be required using differential-equation formulations. In all cases an excellent agreement with the expected results from theory was obtained, confirming the validity and high accuracy of the formulation. On the other hand, there were no a priori conditions included into the formulation or the solution of these problems. On the contrary, results were naturally obtained fulfilling Maxwell’s equations. Indeed, they constitute a further verification of the extraordinary properties of LHM.
This work was supported by the Spanish Government (projects TEC2008-06714-C02-01, TEC2008-06714-C02-02, and CONSOLIDER-INGENIO2010 CSD2008-00068), and by Xunta de Galicia (project IN-CITE08PXIB322250PR).
References and links
1. V. G. Veselago, “The electrodynamics of substances with simultaneously negative values of ε and μ,” Soviet Physics USPEKI 10(4), 509–514 (1968). [CrossRef]
3. N. Engheta and R. W. Ziolkowski, “A positive future for double-negative metamaterials,” IEEE Trans. Microw. Theory Tech. 53(4), 1535–1556 (2005). [CrossRef]
4. C. D. Moss, T. M. Grzegorczyk, Y. Zhang, and J. A. Kong, “Numerical studies of left-handed materials,” Prog. Electromagn. Res. 35, 315–334 (2002). [CrossRef]
5. P. F. Loschialpo, D. W. Forester, D. L. Smith, F. J. Rachford, and C. Monzon, “Optical properties of an ideal homogeneous causal left-handed material slab,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 70(3), 036605 (2004). [CrossRef] [PubMed]
6. Z. G. Dong, S. N. Zhu, H. Liu, J. Zhu, and W. Cao, “Numerical simulations of negative-index refraction in wedge-shaped metamaterials,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 72(1), 016607 (2005). [CrossRef] [PubMed]
7. R. S. Schechtera and S. T. Chun, “Large finite-difference time domain simulations of a left-handed metamaterial lens with wires and resonators,” Appl. Phys. Lett. 91(15), 154102 (2007). [CrossRef]
8. Ö. Ergül, T. Malas, C. Yavuz, A. Ünal, and L. Gürel, “Computational Analysis of complicated metamaterial structures using MLFMA and nested preconditioners,” in Proceedings of European Conference on Antennas and Propagation (EuCAP), Edinburgh, UK, 11–16 Nov. 2007.
9. F. Olyslager, L. Meert, and K. Cools, “The fast multipole method in electromagnetics applied to the simulation of metamateriales,” J. Comput. Appl. Math. 215(2), 528–537 (2008). [CrossRef]
10. P. Yla-Oijala, Ö. Ergül, L. Gürel, and M. Taskinen, “Efficient surface integral equation methods for the analysis of complex metamaterial structures,” in Proceedings of European Conference on Antennas and Propagation (EuCAP), Berlin, Germany, 23–27 Mar. 2009.
11. L. Gürel, Ö. Ergül, A. Ünal, and T. Malas, “Fast and accurate analysis of metamaterial structures using the multilevel fast multipole algorithm,” Prog. Electromagn. Res. 95, 179–198 (2009). [CrossRef]
12. M. Lapine and S. Tretyakov, “Contemporary notes on metamaterials,” IET Microw. Ant. Propag. 1, 3–11 (2007). [CrossRef]
13. D. R. Smith, S. Schultz, P. Markos, and C. M. Soukoulis, “Determination of effective permittivity and permeability of metamaterials from reflection and transmission coefficients,” Phys. Rev. B 65(19), 195104 (2002). [CrossRef]
14. X. Chen, T. M. Grzegorczyk, B. I. Wu, J. Pacheco, and J. A. Kong, “Robust method to retrieve the constitutive effective parameters of metamaterials,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 70(1), 016608 (2004). [CrossRef] [PubMed]
15. D. R. Smith and J. B. Pendry, “Homogenization of metamaterials by field averaging,” J. Opt. Soc. Am. B 23(3), 391–403 (2006). [CrossRef]
16. C. Caloz, C. C. Chang, and T. Itoh, “Full-wave verification of the fundamental properties of left-handed materials in waveguide configurations,” J. Appl. Phys. 90(11), 5483 (2001). [CrossRef]
19. P. F. Loschialpo, D. L. Smith, D. W. Forester, and F. J. Rachford, “Electromagnetic waves focused by a negative-index planar lens,” Phys. Rev. E, Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top. 67, 025602 (2003).
20. C. G. Parazzoli, R. B. Greegor, J. A. Nielsen, M. A. Thompson, K. Li, A. M. Vetter, M. H. Tanielian, and D. C. Vier, “Performance of a negative index of refraction lens,” Appl. Phys. Lett. 84(17), 3232–3234 (2004). [CrossRef]
22. B. J. Justice, J. J. Mock, L. Guo, A. Degiron, D. Schurig, and D. R. Smith, “Spatial mapping of the internal and external electromagnetic fields of negative index metamaterials,” Opt. Express 14(19), 8694–8705 (2006). [CrossRef] [PubMed]
23. R. F. Harrington, Field Computation by Moment Method. (NY: IEEE Press, 1993).
24. L. Landesa, J. M. Taboada, J. L. Rodriguez, F. Obelleiro, J. M. Bertolo, J. C. Mouriño, and A. Gomez, “Analysis of 0.5 Billion Unknowns Using a Parallel FMM-FFT Solver,” in Proceedings of IEEE Antennas and Propagation Society International Symposium, Charleston, SC, USA, 1–5 Jun. 2009.
25. J. M. Taboada, L. Landesa, F. Obelleiro, J. L. Rodriguez, J. M. Bertolo, M. G. Araujo, J. C. Mouriño, and A. Gomez, “High scalability FMM-FFT electromagnetic solver for supercomputer systems,” IEEE Ant. Propag. Mag. 51(6), 20–28 (2009). [CrossRef]
26. M. G. Araújo, J. M. Taboada, F. Obelleiro, J. M. Bértolo, L. Landesa, J. Rivero, and J. L. Rodríguez, “Supercomputer aware approach for the solution of challenging electromagnetic problems,” Prog. Electromagn. Res. 101, 241–256 (2010). [CrossRef]
27. A. J. Poggio, and E. K. Miller, Computer Techniques for Electromagnetics (Elmsford, NY: Permagon, 1973), Chap. 4.
28. S. M. Rao and D. R. Wilton, “E-field, H-field, and combined field solution for arbitrarily shaped three-dimensional dielectric bodies,” Electromag. 10(4), 407–421 (1990). [CrossRef]
29. P. Ylä-Oijala and M. Taskinen, “Application of combined field integral equation for electromagnetic scattering by dielectric and composite objects,” IEEE Trans. Antenn. Propag. 53(3), 1168–1173 (2005). [CrossRef]
30. P. Ylä-Oijala, M. Taskinen, and S. Järvenpää, “Surface integral equation formulations for solving electromagnetic scattering problems with iterative methods,” Radio Sci. 40(6), RS6002 (2005). [CrossRef]
31. P. Ylä-Oijala, M. Taskinen, and J. Sarvas, “Surface integral equation method for general integral equation method for general composite metallic and dielectric structures with junctions,” Prog. Electromagn. Res. 52, 81–108 (2005). [CrossRef]
32. Ö. Ergül and L. Gürel, “Comparison of integral-equation formulations for the fast and accurate solution of scattering problems involving dielectric objects with the multilevel fast multipole algorithm,” IEEE Trans. Antenn. Propag. 57(1), 176–187 (2009). [CrossRef]
33. S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” IEEE Trans. Antenn. Propag. 30(3), 409–418 (1982). [CrossRef]
34. D. R. Wilton, S. M. Rao, A. W. Glisson, D. H. Schaubert, O. M. Al-Bundak, and C. M. Butler, “Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains,” IEEE Trans. Antenn. Propag. 32(3), 276–281 (1984). [CrossRef]
35. R. E. Hodges and Y. Rahmat-Samii, “The evaluation of MFIE integrals with the use of vector triangle basis functions,” Microw. Opt. Technol. Lett. 14(1), 9–14 (1997). [CrossRef]
36. R. D. Graglia, “On the numerical integration of the linear shape functions times the 3-D green’s function or its gradient on a plane triangle,” IEEE Trans. Antenn. Propag. 41(10), 1448–1455 (1993). [CrossRef]
37. P. Yla-Oijala and M. Taskinen; “Calculation of CFIE impedance matrix elements with RWG and, ” IEEE Trans. Antenn. Propag. 51(8), 1837–1846 (2003). [CrossRef]