Open Access
Add to BibSonomyAbstract
Achieving negative permittivity and negative permeability signifies a key topic of research in the design of metamaterials. This paper introduces a level-set based topology optimization method, in which the interface between the vacuum and metal phases is implicitly expressed by the zero-level contour of a higher dimensional level-set function. Following a sensitivity analysis, the optimization maximizes the objective based on the normal direction of the level-set function and induced current flow, thereby generating the desirable patterns of current flow on metal surface. As a benchmark example, the U-shaped structure and its variations are obtained from the level-set topology optimization. Numerical examples demonstrate that both negative permittivity and negative permeability can be attained.
©2010 Optical Society of America
1. Introduction
Metamaterials are a new class of man-made materials possessing outstanding physical properties unattainable from any known natural substance, which have enjoyed a significant success in diverse areas of research recently. Those supernatural properties (e.g. negative refractive index) make some extraordinary applications, such as invisible cloak [1] and super-lens [2], physically possible. In this respect, Russian physicist Veselago made an enlightening contribution [3] to modern metamaterials as early as 1960s, while a strong interest in metamaterials was not sparked until the groundbreaking studies by Pendry and Smith in late 1990s [2, 4, 5].
Metamaterials are usually constructed by periodically-repeated base cells, which are used as a representative volume element for in-depth studies. It is recognized that three key design factors determine the properties of metamaterials, namely the size, composition and configuration of base cell. Of these three factors, the configuration appears the most critical as metamaterials gain their properties from electromagnetic resonance, which can be triggered only in some special microstructures. The most commonly-adopted shape for base cell is the split-ring resonator (SRR) [4], in which the resonance takes place in the structural gaps. On the basis of SRR, more sophisticated configurations like double SRRs, L/U-shaped structures have been devised to construct different metamaterials. Recently, the genetic algorithm, one of the typical structural optimization methods [6], has been used to tailor the optical properties of metamaterials in the negative territory [7–9]. These specially-designed base-cell configurations indeed result in the desired negative permittivity and permeability of metamaterials in certain radiation frequencies [10, 11].
This study concerns the configuration design of metamaterials, which can be considered as a problem of distributing the given basis metal material within the base cell (namely the design domain), such that some specific objectives can be achieved under certain constraints. Structural topology optimization [6] is an effective design tool for finding the optimal configurations for a full range of engineering problems [6, 12–14], which has been well developed to address this issue. In the traditional density-based topology optimization, the basis material is allocated in an element-wise fashion within a finite element framework, in which the local relative density (1/0) is optimized for an ideal state (i.e. 1–solid or 0–void). Since the alternative 1/0 density needs to be relaxed for variable continuity in some topological optimization algorithms, intermediate elemental densities leading to an ambiguous representation appear unavoidable in a final design [15]. Furthermore, the elemental density representation inevitably results in zig-zag boundaries that are unfavorable to some boundary-sensitive structures, e.g. fluidic and electromagnetic problems, which naturally require smoothly-curved boundaries. To tackle such problems, the level-set method [16] was introduced to model the structure through the continuous zero-level contour of a higher dimensional level-set function. It allows smooth tracking of dynamically-changing interfaces, regardless of how fast or sophisticated the changes could be. The level-set technique has been successfully applied to many engineering problems like structural topology optimization [17–19] and inverse scattering problems (shape reconstruction from scattered waves) [20, 21]. However, it is noted that limited studies are available on devising metamaterials using topology optimization and none of them has utilized the level-set method.
This paper aims to introduce the level-set method into metamaterials design by effectively controlling the induced current flow on a metal surface. To simplify the optimization, the governing equation used for calculating the current flow is the electric field integral equation (EFIE) [22] rather than the typical Maxwell's equations. The normal direction of the level-set function is used to define the sign of the current intensity. The integral of the sign-weighted values over the metal surface is maximized throughout the optimization to obtain a certain circuit current pattern that is required in the typical SRR design [4] and other novel metamaterials at their resonant frequencies [23].
To drive the evolution of vacuum/metal interface towards an optimum, the sensitivity of the objective function with respect to the variation of metal-vacuum interface is derived from the shape derivative and an adjoint variable method. This sensitivity is subsequently used as the normal velocity in the level-set procedure. Unlike conventional size or shape optimization adopted in metamaterial design, the level-set method presented herein is capable of optimizing the size, shape and topology of the base cell altogether. As the benchmark examples, several U-shaped structures are reproduced from the level-set optimization framework proposed. The numerical experiments demonstrate that the optimal configurations do attain negative permittivity and negative permeability simultaneously.
2. Materials and methods
For a 2D metal structure subjected to an incident electric field , the induced surface current flow is governed by EFIE [22] as
where Ω stands for the closed region occupied by the metal material in the base cell. The magnetic vector potential and scalar potential V are defined as and , respectively, with the surface charge density equal to and Green's function , where , is the surface divergence and defines the distance between an observation point and a source point . The wave number is , where ω and λ 0 refer to the frequency and wavelength of the incident wave, and ε 0 and μ 0 denote the permittivity and permeability of vacuum. It must be pointed out that, in this paper, we use EFIE rather than the Maxwell's equations to simplify the optimization process. Such a simplification is sensible because EFIE can be derived directly from the Maxwell's equations.Unlike conventional density-based topology optimization [15, 24], the level-set method implicitly represents the metal configuration in the negative territory of a higher-dimensional level-set function, namely [16]. For example, the green (light grey in black-white print) region in Fig. 1(a) is the zero-level slice of the level-set function depicted in Fig. 1(b). Figure 1 illustrates that some typical metamaterial configurations (e.g. double SRR here) can be mathematically represented via such a level-set model.
The optimal evolution of the metal phase can be traced mathematically via the well-known Hamilton–Jacobi (HJ) equation,
as long as the normal velocity (i.e. the gradient of the objective function) can be determined.The resonant response to the incident electrical field is the main reason why metamaterials have extraordinary electromagnetic properties [1, 2, 4, 5, 23, 25]. Pendry et al. [4] showed that the magnetic field induced by the surface current flow changes the original magnetic field and results in negative effective permeability of a cylindrical wire. It implies that the induced current flow plays a key role in determining the electromagnetic properties for metamaterials. Figures 2 (b1)-(b3) illustrate the surface current flow (obtained by solving the Maxwell's equations) in several typical configurations (Figs. 2 (a1)-(a3)), where the induced electric field points upward. From these sub-figures, it can be seen that the effective permittivity (Figs. 2(c1)-(c3)) and permeability (Figs. 2 (d1)-(d3)) can attain negative values when a current circuit is generated on the metal surface (e.g. the counterclockwise flow for the well-known SRR structure in Fig. 2(b2)). It is noted that the passive metamaterial can always have positive real part of impedance and positive imaginary part of refractive index; while their real part of permittivity and imaginary part of permeability can attain the negative values in a certain frequency range as per their relation to the impedance and refractive index. For this reason, we only plot these two more indicative parts (i.e. Re(ε), and Im(μ)) herein. Based on this, Liu et al. [23] observed a circular current flow in a U-shaped structure and inferred that such a pattern of current flow could lead to two resonances, which were dominantly associated with plasmon excitations in its two arms. It is noted that the effective properties for these and following structures are obtained from a so-called S-parameter retrieval method [26] via printing the metal material (gold, with 36 μm outer dimension in a lattice of 50 μm) on a 4μm thick substrate (semi-insulating gallium arsenide GaAs). Such a modelling technique has been verified in literature, experimentally and numerically [23].
Fig. 2 The effective electromagnetic properties of different structures and their induced surface current flow: (a1)-(a3) structural configurations; (b1)-(b3) induced surface current flow; (c1)-(c3) effective permittivity; (d1)-(d3) effective permeability.
It is noted that the metal surface is constrained to a square with 36 μm width centrally located in the periodically-repeated base cell of a square with 50 μm width, which means that there is a 14 μm gap between the metal structures in two neighboring base cells. In other words, the metal region actually does not touch the boundary of base cell, and thus there is no need to prescribe periodicity in the metal edge when solving EFIE. However, in the level-set model that defines base cell, periodic boundary conditions must be applied. For this purpose, we applied additional ghost nodes surrounded the base cell, as did in [19].
On the other hand, a symmetric current flow may not necessarily result in the desired negative permittivity and permeability. For instance, if we rotate the electric field by 90° with respect to the split in SRR, as shown in Fig. 2(a3), the permittivity barely attains negative territory while the permeability keeps positive (Figs. 2(c3)-(d3)). In this situation, the desired circuit current does not take place. A similar symmetric current flow with regard to the incident electric field (Fig. 2(b1)) also appears on the surface of the circular ring structure (Fig. 2(a1)), which leads to positive permittivity and permeability at all frequencies (Figs. 2(c1)-(d1)). In fact, Linden et al. [25] experimentally proved that a symmetric current flow pattern is unable to render negative permittivity/permeability. The effective properties in Fig. 2 are extracted in line with the method proposed in [27] by solving the Maxwell's equations, which approximating the periodic boundaries by the setting the perfect electric and perfect magnetic boundaries on the two sets of opposite boundaries parallel to the propagating direction.
The abovementioned studies provide strong evidence that the circuit current flow plays a critical role in determining the characteristics of metamaterials. To achieve such a favorable pattern of current flow, we thus formulate an objective function as
where the signed weighting factor when the local current flow is counterclockwise ( + 1) or clockwise (–1), respectively. The unit vector in S is obtained from the projection of the normal direction of the 3D level-set function to the 2D metal domain Ω. This objective function tends to be zero for a nearly-symmetric pattern of current flow, such as those in Figs. 2 (b1) and (b3), because the signed weighting factor S takes –1 and + 1, respectively, on the two vertical sides; thus the current intensity (the norm of current flow) is largely cancelled out. While for the circuit current flow in Fig. 2 (b2), this objective becomes highly positive as the direction of current flow keeps constant with respect to the normal direction (i.e. S keeps + 1 all the time). The cost function (in both magnitude and sign) plays a key role on driving the topology optimization in terms of velocity (sensitivity) in the proposed level-set procedure. It is also a measure of the optimization extent. In fact, even though the sign of the current flow is desirably right, the metamaterial structure can continue evolving toward an optimum. The following examples showed that the magnitude of the current flow can increase hundreds of times throughout the optimization process. In some parts of the optimal structure, the current flow can attain rather great values.The EFIE system is solved by the method of moment (MoM) [28], in which the metal region is adaptively discretized into triangular elements (Fig. 3(a) ) in terms of the local value of the level-set function. Such an adaptive triangular mesh can capture the curved boundaries rather precisely with relatively low computing cost. The points on the zero-level contour of the level-set function are selected as the triangular vertices in the mesh thus zig-zag boundaries are avoided. After solving the EFIE system, the current values at the central nodes of the triangular elements are interpolated to the nodes of the rectangular mesh that will be used in the level-set algorithm. While the rectangular mesh does not changes after it is generated at the first step as its covering region, namely the base cell, is fixed in the optimization. However, the triangular mesh should be regenerated in line with the topological variation of the metal surface. Following this, its according level-set values need to be remapped in each step.
Fig. 3 The current flow obtained by solving EFIE: (a) the adaptive triangular mesh for a SRR; (b) the current flow on a SRR; (c) the current flow on a circular ring without aperture.
The current flows obtained from EFIE, Figs. 3 (b)-(c), exhibit the similar distributions to those obtained from the Maxwell's equations (Figs. 2 (b2) and (b1), respectively). Such equivalence provides a justification of replacing the Maxwell's equations with EFIE in this paper.
The derivative of the objective function with respect to the metal shape Ω is obtained from the shape derivative and the adjoint variable method, given by
with and the adjoint variable is the solution to the adjoint equation of EFIE, given in the weak form bywhere denotes a virtual variable.In a level-set model, its normal velocity can be expressed as
Along with this direction, the shape of the base cell can evolve towards a state, enabling the objective function to increase, mathematically given by3. Results and discussion
Three examples, each starting from different initial configurations, are presented herein to demonstrate how the level-set procedure optimizes the shape and topology of base-cell structures. The base-cell models in the following examples are subjected to an incident wave with unit vertical electrical field. As needed in metamaterials, the wavelength of the incident electric field can be 10 times larger than the size of the base cell. The level-set function should be reinitialized in each iteration to keep itself as a sign-distanced function. In this paper, the rectangular mesh used for the level-set algorithm is 100×100, while the triangular mesh for MoM is adaptive. To make the structures symmetrical as commonly-used in metamaterials, we imposed symmetrical constraint to the base cell structures during the optimization. Nevertheless, our numerical tests showed that this method can achieve nearly-symmetrical structures even without this symmetrical constraint.
Example 1 starts from a single circle located in the centre of the base-cell domain. Figure 4 shows that the objective value raises slowly when the metal region evolves from a full circle (Fig. 4(a)) to a semi-circle (Fig. 4(b)) over the first 30 iterations (the red dots on the curve correspond to the base-cell structures on the top of Fig. 4, from (a) to (e)). Following the appearance of a dent on the left edge in iteration 60 (Fig. 4(c)), the objective value rises gradually as this dent becomes deeper (Fig. 4(d)). Subsequently, the objective value rises sharply to a level much higher than that of the initial value and the design finally evolves into a U-shaped structure (Fig. 4(e)).
It must be pointed that the cost functional could increase infinitely provided that the discretization of the local mesh could be infinitely small. In other words, the bottom of the U-shape structure should be as thin as possible to make the induced current as high as possible, theoretically. In practice, however, the minimum size of geometric feature is limited by the fabrication technologies available. Besides, the thinner the metal part, the easier to be burnt out by the current flow. Thus, the mesh density in such critical region with most interest cannot be increased further. With a limited minimum mesh size, the structure will break at the part with maximal current flow, and accordingly, the objective value will plunge beyond Fig. 4(e), leaving another round of optimization starting from these separated parts. Thus, such factors to a certain extent restrict the maximum value that the objective functional can achieve.
In Example 1, the real part of the permittivity (Re(ε)) and imaginary part of permeability (Im(μ)) are plotted in Figs. 5(a)-(b) , respectively, illustrating that such a U-shaped structure has both decent negative Re(ε) and Im(μ) at around 0.9 Tera Hertz (THz). On the contrary, the Re(ε) of the initial design (circular shape) attains only a shallow negative territory when the incident frequency f>1.2 THz, while its Im(μ) remains invariably positive in the entire frequency range considered.
Fig. 5 The effective electromagnetic properties for structures in initial and final stages: (a) the real part of permittivity (Re(ε)); (b) the imaginary part of permeability (Im(μ)). (Example 1)
Although a significant geometric change is resulted, this example is essentially subjected to a shape optimization without topological change occurring during the design process. It must be noted that as a benchmark example, such a U-shaped structure has been well studied in metamaterial design. It has been fabricated in 2D and 3D, and does have negative electromagnetic properties even at visible frequencies [10, 29].
In Example 2, the design starts from two identical circles (Fig. 6(a) , named structure 0 at the 0th iteration step). The optimization makes these two circles merge in an early stage, shaping two dents on bilateral edges (Fig. 6(b), structure 5). As the design progresses, the right dent gradually disappears, while the left dent deepens further (Fig. 6(c), structure 50). Finally, the growth of the left dent also shapes the metal region into a U-shaped structure (Fig. 6(d), structure 75), agreeing well with the previous benchmarking example.
Fig. 6 The metal shapes in different stages for Example 2 starting from two identical circles: (a) structure 0; (b) structure 5; (c) structure 50; (d) structure 75.
The effective electromagnetic properties for these several snapshots in Fig. 6 are plotted in Fig. 7 , where both Re(ε) and Im(μ) attain negative territory at around 1.1 THz for the optimized U-shaped structure. Figure 7 also exhibits that the negative peaks of Re(ε) deepen for the corresponding structures from Figs. 6(a) to (d), while the positive peaks for Re(ε) and Im(μ) drift toward the left during the optimization. In this example, topological change takes place (the mergence of two circles). It seems that the design progresses much faster in the topology optimization compared with the pure shape optimization given in the previous example (only 75 iterations herein compared with 120 iterations in Example 1).
Fig. 7 The effective electromagnetic properties for structures in different stages: (a) the real part of permittivity (Re(ε)); (b) the imaginary part of permeability (Im(μ)). (Example 2)
Example 3 starts from an initial structure that consists of five circular holes in a big disk (Fig. 8(a) ). It is observed that due to its topological sophistication, more frequent topology changes occur during the optimization, accompanying the mergence and isolation of the holes. With two sets of controlling parameters, two different final structures, Figs. 8(b)-(c), are obtained after 32 and 40 iterations, respectively.
Fig. 8 The initial structure (a) and different final structures (b-c) with different control parameters (Example 3).
Interestingly, Fig. 9 indicates that these two optimal configurations have very similar effective electromagnetic properties even though their final shapes look somewhat different. It is observed that the effect of the counterclockwise current flow on the major right metal region is weakened by the clockwise flow on the isolated small left object in Figs. 8(b)-(c). This is why the minimal Re(ε) and Im(μ) are shallower than those in Examples 1-2. Similarly, the depth and width of the U-shape play a key role on determining the effective values. These factors can be used to tune the electromagnetic properties accordingly.
Fig. 9 The effective electromagnetic properties for the structures in different design stages: (a) the real part of permittivity (Re(ε)); (b) the imaginary part of permeability (Im(μ)). (Example 3)
4. Conclusions
To devise the base-cell configuration for metamaterials subject to an incident electrical field, we introduced a level-set based topology optimization algorithm. In this method, a scalar level-set function is used to implicitly express the metal configuration via its zero-level contour. The objective function is defined in terms of the current norm and the normal direction of the level-set function. The shape derivative and adjoint variable method are adopted to derive the sensitivity of the objective function with respect to shape variation. From the level-set procedure presented, several different U-shaped structures are obtained. Their permittivity and permeability attain negative values within certain frequency ranges. It must be pointed out that the level-set framework presented herein allows accommodating other objective functions relevant to metamaterials design.
Acknowledgments
This study has been supported by Australian Research Council (ARC). We thank Joseph Cadman at The University of Sydney for his comments and a proof reading of the manuscript.
References and links
1. J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science 312(5781), 1780–1782 (2006). [CrossRef] [PubMed]
2. J. B. Pendry, “Negative refraction makes a perfect lens,” Phys. Rev. Lett. 85(18), 3966–3969 (2000). [CrossRef] [PubMed]
3. V. G. Veselago, “The electrodynamics of substances with simultaneously negative values of ε and μ,” Sov. Phys. USPEKI 10(4), 509–514 (1968). [CrossRef]
4. J. B. Pendry, A. J. Holden, D. J. Robbins, and W. J. Stewart, “Magnetism from conductors and enhanced nonlinear phenomena,” IEEE Trans. Microw. Theory Tech. 47(11), 2075–2084 (1999). [CrossRef]
5. D. R. Smith, W. J. Padilla, D. C. Vier, S. C. Nemat-Nasser, and S. Schultz, “Composite medium with simultaneously negative permeability and permittivity,” Phys. Rev. Lett. 84(18), 4184–4187 (2000). [CrossRef] [PubMed]
6. M. P. Bendsøe, and O. Sigmund, Topology Optimisation: theory, methods, and applications (Springer, Berlin; New York, 2003).
7. D. H. Kwon and D. H. Werner, “Low-index metamaterial designs in the visible spectrum,” Opt. Express 15(15), 9267–9272 (2007). [CrossRef] [PubMed]
8. P. Y. Chen, C. H. Chen, H. Wang, J. H. Tsai, and W. X. Ni, “Synthesis design of artificial magnetic metamaterials using a genetic algorithm,” Opt. Express 16(17), 12806–12818 (2008). [CrossRef] [PubMed]
9. J. A. Bossard, S. Yun, D. H. Werner, and T. S. Mayer, “Synthesizing low loss negative index metamaterial stacks for the mid-infrared using genetic algorithms,” Opt. Express 2009(17), 14771–14779 (2009). [CrossRef]
10. N. Liu, H. C. Guo, L. W. Fu, S. Kaiser, H. Schweizer, and H. Giessen, “Three-dimensional photonic metamaterials at optical frequencies,” Nat. Mater. 7(1), 31–37 (2008). [CrossRef]
11. V. M. Shalaev, W. Cai, U. K. Chettiar, H. K. Yuan, A. K. Sarychev, V. P. Drachev, and A. V. Kildishev, “Negative index of refraction in optical metamaterials,” Opt. Lett. 30(24), 3356–3358 (2005). [CrossRef]
12. Q. Li, G. P. Steven, O. M. Querin, and Y. M. Xie, “Shape and topology design for heat conduction by evolutionary structural optimisation,” Int. J. Heat Mass Transfer 42(17), 3361–3371 (1999). [CrossRef]
13. G. P. Steven, Q. Li, and Y. M. Xie, “Evolutionary topology and shape design for general physical field problems,” Comput. Mech. 26(2), 129–139 (2000). [CrossRef]
14. S. W. Zhou and Q. Li, “The relation of constant mean curvature surfaces to multiphase composites with extremal thermal conductivity,” J. Phys. D Appl. Phys. 40(19), 6083–6093 (2007). [CrossRef]
15. M. Zhou and G. I. N. Rozvany, “The COC algorithm. II: Topological, geometrical and generalized shape optimization,” Comput. Methods Appl. Mech. Eng. 89(1-3), 309–336 (1991). [CrossRef]
16. S. Osher and J. A. Sethian, “Front propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations,” J. Comput. Phys. 79(1), 12–49 (1988). [CrossRef]
17. G. Allaire, F. Jouve, and A. M. Toader, “Structural optimization using sensitivity analysis and a level-set method,” J. Comput. Phys. 194(1), 363–393 (2004). [CrossRef]
18. M. Y. Wang, X. M. Wang, and D. M. Guo, “A level set method for structural topology optimization,” Comput. Methods Appl. Mech. Eng. 192(1-2), 227–246 (2003). [CrossRef]
19. S. W. Zhou and Q. Li, “A variational level set method for the topology optimization of steady-state Navier-Stokes flow,” J. Comput. Phys. 227(24), 10178–10195 (2008). [CrossRef]
20. M. Burger and S. J. Osher, “A survey on level set methods for inverse problems and optimal design,” Eur. J. Appl. Math. 16(2), 263–301 (2005). [CrossRef]
21. O. Dorn and D. Lesselier, “Level set methods for inverse scattering,” Inverse Probl. 22(4), R67–R131 (2006). [CrossRef]
22. A. W. Maue, “On the formulation of a general scattering problem by means of an integral equation,” Z. Phys. 126, 601–618 (1949). [CrossRef]
23. W. J. Padilla, M. T. Aronsson, C. Highstrete, M. Lee, A. J. Taylor, and R. D. Averitt, “Electrically resonant terahertz metamaterials: Theoretical and experimental investigations,” Phys. Rev. B 75, 041102 041101/041104 (2007). [CrossRef]
24. Y. M. Xie and G. P. Steven, “A simple evolutionary procedure for structural optimization,” Comput. Struc. 49(5), 885–896 (1993). [CrossRef]
25. S. Linden, C. Enkrich, M. Wegener, J. Zhou, T. Koschny, and C. M. Soukoulis, “Magnetic response of metamaterials at 100 Terahertz,” Science 306(5700), 1351–1353 (2004). [CrossRef] [PubMed]
26. D. R. Smith, D. C. Vier, T. Koschny, and C. M. Soukoulis, “Electromagnetic parameter retrieval from inhomogeneous metamaterials,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 71(33 Pt 2B), 036617 (2005). [CrossRef] [PubMed]
27. G. Lubkowski, R. Schuhmann, and T. Weiland, “Extraction of effective metamaterial parameters by parameter fitting of dispersive models,” Microwave Optical Tech. Lett. 49(2), 285–288 (2007). [CrossRef]
28. R. F. Harrington, Field Computation by the Moment Methods (IEEE Press, New York, 1993).
29. C. Enkrich, M. Wegener, S. Linden, S. Burger, L. Zschiedrich, F. Schmidt, J. F. Zhou, T. Koschny, and C. M. Soukoulis, “Magnetic metamaterials at telecommunication and visible frequencies,” Phys. Rev. Lett. 95(20), 203901 (2005). [CrossRef] [PubMed]
