## Abstract

We reformulate the problem of a steerable left-handed antenna reported by Matsuzawa et al. [IEICE Trans. Electron. **E89-C**, 1337 (2006)] from the view point of structural electromagnetic resonance of the unit structure. We show that there are two such resonances with different spatial symmetries in the relevant frequency range, which result in the formation of two electromagnetic bands with opposite signs of curvature at the Γ point of the Brillouin zone. We derive an expression of dispersion curves based on the tight-binding picture and show that the dispersion of the two bands is linear in the vicinity of the Γ point in the case of accidental degeneracy only if the symmetry of the two resonance states satisfies certain conditions. We also show that the refraction angle can be designed by changing the lattice constant of the arrayed unit structures, since the band width is modified due to the change in the electromagnetic *transfer integral*.

© 2010 Optical Society of America

## 1. Introduction

Negative refraction realized by left-handed materials, or metamaterials, has attracted a great deal of both scientific and technological attention [1–7]. As a particular application, a microwave steerable antenna based on negative refraction was proposed and demonstrated [8–10]. Microwave antennas of this kind are important, since they may be used for automotive radar sensors for adaptive cruise control and pre-crash safety systems [11]. Their previous analysis and design as left-handed materials were successful and described the properties of the steerable antenna well.

However, from a physical point of view, the relation between the physical substance and the antenna properties is not necessarily clear. Especially unclear is the physical origin of the two electromagnetic bands with opposite signs of curvature at the Γ point of the Brillouin zone, which play the most essential role in the steerability. Solving these problems is important not only for practical applications but also for fundamental physics, since the connection between the behavior of the electromagnetic field in the sub-wavelength scale and the electromagnetic response of the metamaterial as a whole is an essential issue.

Generally speaking, metamaterials utilize local electromagnetic resonance to realize their unique properties such as negative refractive index. In the microwave region, structural resonances are used for this purpose, whereas plasmonic resonance is also used in the optical region. Because the resonant state is usually well localized on the metallic unit structure, the tight-binding approximation may give an accurate description of the electromagnetic field of the total system.

If the metamaterial consists of a regular array of unit structures as in the case of the steerable antenna analyzed in Ref. [8], the electromagnetic field yields a Bloch band. This situation is similar to the case of photonic crystals [12,13], although the specimens are usually composed of dielectric materials with a positive dielectric constant in the latter case so that local confinement of the electromagnetic field is weak and the Bloch modes can be well described by linear combination of extended plane waves. In contrast, the electromagnetic Bloch modes in the former case are far from plane waves.

The Bloch bands of metallic photonic crystals have also been analyzed by several authors. The simplest case may be a regular two-dimensional array of metallic circular cylinders [14], in which the Mie resonance states play the role of localized states that give the basis functions for the tight-binding approximation. However, the simple tight-binding approximation that is used for electronic band calculation is not applicable, since the amplitude of resonant states follows the power-law dependence on the distance from the scatterer so that their wave functions are not normalizable. To resolve this problem, we may impose a periodic boundary condition on the wave functions of the resonant state in a sufficiently large volume compared with the size of the scatterer and normalize them in this volume as we usually do for plane waves. In Ref. [14], the Bloch bands were actually calculated by the finite-difference time-domain (FDTD) method [15, 16] assuming the Drude dielectric function, whereas their symmetries were predicted by group theory based on the tight-binding picture and they agreed with the results of the FDTD calculation well.

Thus in this paper, we reformulate the problem of the steerable antenna without the idea of a negative refractive index. Rather, we simply calculate the structural resonances of the unit structure by the FDTD method and investigate their consequences first qualitatively by the tight-binding picture and later by quantitative numerical calculation. In Section 2, based on the tight-binding picture, we show the presence of two Bloch bands with opposite signs of curvature at the Γ point by assuming two resonant states localized on a unit structure with different parities, and prove the steering property. In Section 3, we first show by the FDTD method that the unit structure of the specimen actually has two structural electromagnetic resonances in the relevant frequency range whose symmetries are mutually different. Next, we calculate the dispersion relation of eigenmodes to show that the two relevant bands have opposite signs of curvature at the Γ point. In addition, we show by the tight-binding calculation and group theory that the two dispersion curves have finite slopes at the Γ point in the case of accidental degeneracy if their symmetries satisfy a certain condition. Finally, we present a brief discussion on designing the steerable antenna by adjusting device parameters. A summary is given in Section 4.

## 2. Theory

Figure 1 shows the structure of the specimen analyzed in Ref. [8] and this paper. The specimen consists of a one-dimensional regular array of metallic unit structures fabricated on a dielectric slab with a ground electrode on the back surface. Both the unit structure and the regular array as a whole have the *C*_{2v} symmetry. Since the antenna was designed to operate at 76 GHz (free-space wavelength = 3.95 mm), the size of the unit structure is of the order of 1 mm and the period of the array, or the lattice constant *a*, is 0.6 mm. The antenna structure is periodic in the *x* direction. An incident wave with angular frequency *ω _{i}* is supplied in the positive

*x*direction to propagate in the dielectric slab.

Now, let us examine the functionality of the steerable antenna qualitatively by assuming two electromagnetic resonances for the single unit structure whose electric fields have opposite parities with respect to the *x* coordinate. The presence of such resonant states will be proved in Section 3.

Because the resonant states are well localized on the single unit structure, we may apply the tight-binding model to this problem. It is often observed that electronic p states of semiconductors yield energy bands concave-down around the Γ point, while s states yield concave-up energy bands [17]. This difference originates from different parities of p and s states.

The band width *B* in the one-dimensional tight-binding model in the lowest order approximation is given as follows:

*ϕ*(

**r**) is a normalized atomic wave function on the lattice point at the origin,

*ϕ*(

**r**–

**a**) is that on the nearest neighbor lattice point

**a**, and

*U*(

**r**) (< 0) is the attractive potential from all atomic nuclei (ionic cores) except that at the origin [17]. For simplicity, we ignored overlapping of wave functions at the second nearest neighbor and farther lattice points in Eq. (1). The energy band is given by where

*E*

_{0}is the sum of the original electronic energy of an isolated atom and If the wave function is of the s character, that is, if its parity is even, then

*ϕ*(

**r**) and

*ϕ*(

**r**–

**a**) have the same sign at

**r**=

**a**/2 where the dominant contribution is made to the integral of Eq. (1), so that

*B*> 0. On the other hand, if the wave function is of the p character and has the odd parity,

*B*< 0.

By analogy, we may conclude that we have two dispersion curves for the problem of the steerable antenna, one concave-down and the other concave-up, since the parities of the two resonant states are different from each other. Let us denote their dispersion curves as follows:

*f*(

_{l}*f*),

_{u}*A*(

_{l}*A*), and

_{u}*ν*(

_{l}*ν*) are the eigen angular frequency, band width, and central angular frequency of the lower (upper) band, respectively. Here we should note that the square of the eigen frequency takes the role of the electronic energy, since Maxwell’s wave equation is of the second order for the time variable whereas the Schrödinger equation is of the first order. We further assume so that the two bands have the same frequency at the Γ point (

_{u}*k*= 0). This is the second condition that the steerable antenna should satisfy. As will be shown in Section 3, this condition can be realized, for example, by adjusting the lattice constant.

We should note that here we implicitly assumed that the two bands had mutually different symmetries so that they do not repel each other when they come close. If the two resonant states satisfy a certain condition with respect to their symmetry, a more interesting situation is realized: The two bands coincide on the Γ point and repel each other in the rest of the *k* space, which results in the linear dispersion in the vicinity of the Γ point. This case will be described in the next section.

Now, a more realistic expression for the mode frequency can be derived as follows. We denote the position-dependent dielectric constant of the periodic arrayed structure by *ɛ*(**r**), for which we assume the Drude function

*ω*is the angular frequency and

*ω*is the plasma frequency. We assume that magnetic permeability is equal to unity, since we do not deal with magnetic materials. Then, from the Maxwell equations, we have the wave equation for the magnetic field

_{p}**H**(

**r**

*,t*):

*c*is the light velocity in free space. For a single unit structure, we denote its dielectric constant by

*ɛ*(

_{s}**r**) and assume that we know its eigen function

**H**

_{0}(

**r**) with eigen frequency

*ω*

_{0}and that it is normalized:

*V*is the volume on which we impose the periodic boundary condition. Note that

**H**

_{0}is dimensionless by this definition. Then we approximate the eigen function of the periodic system

**H**

*(*

_{k}**r**) by a linear combination of

**H**

_{0}(

**r**–

*n*

**a**) (

*n*= 1, 2, ··· ,

*N*):

**a**is (

*a*, 0, 0). In the lowest order approximation, eigen frequency of the Bloch band

*ω*is given by

_{k}*δχ*(

**r**) is the difference of 1/

*ɛ*(

**r**) between the periodic system and a single unit structure:

*ω*<

*ω*. Performing integration by parts, the electromagnetic

_{p}*transfer integral L*

_{1}is given by

*ɛ*

_{0}is the permittivity of free space. In Eq. (16), we introduced the eigen function of the electric field

**E**

_{0}given by

*S*of volume

*V*, and is vanishing due to the periodic boundary condition. Since a localized mode is a standing wave, its phase is the same everywhere. Thus we can choose

**E**

_{0}(

**r**) as a real function. Because

*δχ*(

**r**) is non-zero only in the region where

*ɛ*(

_{s}**r**) = 1 (air), we finally obtain

Since the lattice constant is sufficiently smaller than the wavelength, the **E**_{0}(**r**) wave function is expected to keep its sign in the region between two adjacent unit structures. Thus we can follow a similar logic as in the case of p and s electronic states in semiconductor. However, we should take into consideration two differences. The first one is the different transformation properties between electronic wave functions and electromagnetic wave functions, or between scalar waves and vector waves. When we denote the mirror reflection of the *x* coordinate about the origin by *σ _{x}*, then the symmetry of an electronic wave function

*ϕ*(

**r**) is expressed as [18]

**E**(

**r**) is given by [18]

We should also note that on the surface of metallic unit structures, the tangential electric field components parallel to the surface are small due to large metallic conductivity. For the present geometry, *E _{z}* is dominant on the metal surface and so is it at the nearest neighbor lattice points. Thus, Eq. (18) can be approximated by

**E**

_{0}has even (odd) parity with respect to the

*x*coordinate,

*L*

_{1}is positive (negative) and the Bloch band is concave-down (up).

Now, let us go back to Eqs. (4) and (5). In this simplest tight-binding picture, the dispersion curves of the system may look like Fig. 2. In addition to the two cosine bands, light lines given by *ω* = ±*ck* are also plotted. Note that the horizontal axis of this figure is the normalized crystalline momentum (divided by *h̄*) in the *x* direction, and is *not* the wave vector in free space.

Now, we examine the situation in which an incident wave with angular frequency *ω _{i}* is propagated in the dielectric slab in the positive

*x*direction, an internal mode with crystalline momentum

*k*is excited, and a diffracted wave is generated. Here we should note that an internal mode with a positive group velocity (=

_{i}*df*/

_{l,u}*dk*), which is equal to the energy velocity [19], can only be excited because of the conservation of energy flow as far as we neglect reflection at the exit end of the device. So, although both positive and negative

*k*’s are present for each

_{i}*ω*, only the positive (negative)

_{i}*k*mode of the upper (lower) band can be excited by the incident wave.

_{i}When we denote the *z* component of the diffracted wave by *k _{z}*,

*k*should be a real number (See Fig. 3). So, diffraction is possible only in the region above the light line as for photonic crystal slabs [13], whose lower and upper bounds are labelled

_{z}*ω*and

_{l}*ω*in Fig. 2. Then, the diffraction angle

_{u}*θ*is given by

^{−1}is from 0 to

*π*.

Frequency dependence of the diffraction angle is shown for some combinations of parameters in Fig. 4. Note that even if the available range of *k* is limited by two light lines (see Fig. 2), the emission angle always spans from 0 to 180 degrees.

In the actual device demonstrated in Ref. [8], the incident frequency *ω _{i}* was fixed. On the other hand, an additional dielectric plate was placed above the metallic array, and the distance between the plate and the array was varied so that the mode frequencies (dispersion curves) were shifted continuously by the perturbation induced by the dielectric plate. Thus the relative position of the incident frequency to the dispersion curves was adjusted and the diffraction angle changed accordingly.

Based on these qualitative considerations of the characteristics of the metamaterial steerable antenna, we performed quantitative numerical calculation by the FDTD method. We wrote paralleled C^{++} codes with the conventional Yee mesh and the convolutional PML (perfectly matched layer) boundary conditions. Since we were dealing with the microwave range, we assumed a low frequency limit of the analytical Drude function for metallic *ɛ*:

*δ*is a positive infinitesimal to assure the causality,

*γ*is the damping constant, and $\sigma ={\varepsilon}_{0}{\omega}_{p}^{2}/\gamma $ is the conductivity. We used 5.8 × 10

^{7}S/m for

*σ*(value for copper).

To discretize the Maxwell equations in the FDTD scheme, we take mutually different temporal mesh points for **E** and **H**, and calculate **H**_{(}_{n}_{+1/2)}_{Δ}* _{t}* and

**E**

_{n}_{Δ}

*, where*

_{t}*n*is an integer and Δ

*t*is the time-discretization unit. When we apply this procedure to Ampere’s law,

Resonance frequencies were obtained by finding peaks in the Fourier transform of the temporal variation of the electromagnetic field after pulsed excitation. The field distributions of the resonant states were calculated by cw (continuous wave) excitation at the resonance frequencies. For the dispersion curves, we imposed the following Bloch boundary conditions on the electric and magnetic fields,

and found the*k*-dependent resonance frequencies by searching peaks in the Fourier transform of the temporal variation of the electromagnetic field after pulsed excitation.

## 3. Results and discussion

Let us examine the resonant states of a single isolated unit structure first. We found two peaks (82.8 GHz and 96.6 GHz) in the Fourier spectrum in the relevant frequency range. Their quality factors were 192 and 59, respectively. Figure 5 shows their field distributions, both of which were localized on the unit structure. It should be noted that the higher-frequency resonance has an odd parity with respect to the *x* coordinate, whereas the lower-frequency resonance has an even parity. As a result, we may assume that the higher band is concave-up and the lower band is concave-down according to the discussion in the previous section, although the assumption of negative *ɛ* does not rigorously apply to the present case. Nonetheless, the above assumption actually happened as shown in Fig. 6, where among the two dispersion curves, the upper one is concave-up and the lower one is concave-down at the Γ point.

These dispersion curves share some features with Fig. 3 of Ref. [8] regarding the number of modes around the Γ point and their frequency ranges. But there are also some differences. First, the slope of the two dispersion curves were flat at the Γ point in our case. This feature was expected from time reversal symmetry, which resulted in the relation

(See, for example, Section 2.3 of Ref. [13]). So, if the band is not degenerate and*ω*is, as usual, a smooth analytic function of

*k*around the Γ point, the Taylor expansion of

*ω*with respect to

*k*has terms of even orders of

*k*only, so that

*dω/dk*= 0 for

*k*= 0.

Second, the two frequencies at the Γ point are not exactly equal to each other. Since the unit structure has the *C*_{2}* _{v}* symmetry, eigenmodes at the Γ point are all one-dimensional, i.e., non-degenerate [18]. So, the coincidence of eigen frequencies of the two modes, if any, is accidental degeneracy, which is brought about by a particular choice of sample parameters. The case of accidental degeneracy will be described later. Thus the probability of having different frequencies is 100 % from the purely mathematical point of view. However, as will also be shown later, we can make them as close to each other as we like by adjusting sample parameters.

Third, we found that a third band was located very close to the light line that had the same parity as the other two bands with respect to the *y* coordinate and showed the anti-crossing behavior that occurred due to intermixing when its frequency crossed the dispersion curves of the two flat bands. The third band originated from the lowest TM (transverse magnetic) mode of the dielectric slab whose dispersion is given by the solution of the following secular equation.

*ɛ*and

_{d}*d*are the dielectric constant and thickness of the slab waveguide, and

*k*and

_{z}*κ*are given by As is also apparent from this anti-crossing behavior, the band modes shown in Fig. 6 can be excited by incident waves of the TM polarization.

Now, let us examine the case of accidental degeneracy. First, we should note that the two resonant states shown in Fig. 5 belong to different irreducible representations of the *C*_{2}* _{v}* point group. By consulting its character table given in Table 1, it is found that the lower- and higher-frequency modes belong to the

*A*

_{1}and

*B*

_{1}representations, respectively. As for the symmetry of the band states, it is the same as the original resonant states on the Γ point (

*k*= 0) and the boundary of the Brillouin zone (

*k*= ±

*π*/

*a*). On the other hand, it only keeps the parity with respect to the

*y*coordinate in the rest of the Brillouin zone [18].

In the present case, we have *A*_{1} and *B*_{1} modes. Thus the symmetries of the two bands on the Γ point are also *A*_{1} and *B*_{1} so that they are mutually different and do not repel each other even when they come close by choosing device parameters. On the other hand, since the two bands have the same even parity with respect to the *y* coordinate for general *k*, they repel each other when they come close. According to the valance of these two factors, linear dispersion curves are realized in the vicinity of the Γ point in the case of accidental degeneracy.

Let us examine this feature by the tight-binding calculation. We denote the lower and higher resonant states by ${\mathbf{H}}_{0}^{(1)}$ and ${\mathbf{H}}_{0}^{(2)}$, respectively. They have the following symmetry:

We assume as before that their wave equations have been solved and denote their eigen angular frequencies by*ω*

_{1}and

*ω*

_{2}:

*V*is the volume on which we imposed the periodic boundary condition. Then the Bloch wave function in the tight-binding picture must be a linear combination of ${\mathbf{H}}_{0}^{(1)}$ and ${\mathbf{H}}_{0}^{(2)}$:

First, let us evaluate the term proportional to 1/*ɛ _{s}* on the left-hand side of Eq. (42). Multiplying
${\mathbf{H}}_{0}^{(1)*}(\mathbf{r})$ and integrating over

*V*, we obtain

*δχ*,

*B*that ${L}_{1}^{(12)}$ and ${L}_{1}^{(21)}$ have opposite signs. Similarly, the right-hand side of Eq. (42) gives $AV{\omega}_{k}^{2}/{c}^{2}N$. By keeping contributions from wave functions on the origin and the nearest neighbor lattice points, we obtain

First, let us consider a special case in which the two bands are degenerate on the Γ point. By setting *k* = 0 in Eq. (51) and assuming the two eigen frequencies are the same, we obtain the following condition for the degeneracy.

*ω*

_{Γ}. Then, the secular equation reduces to

*k*≪

*a*and keep the lowest order terms with respect to

*k*. Thus we obtain

*k*and have finite slopes whose magnitudes are the same and signs are opposite. In particular, the time reversal symmetry described by Eq. (33) holds. These results coincide with previous findings by the theory of transmission lines [10, 20, 21].

On the other hand, in the general non-degenerate case, we obtain from Eq. (51)

Figure 7 shows examples of dispersion curves with some combinations of device parameters. When we increase ${L}_{0}^{(11)}$ from (a) to (d), the two bands touch each other on the Γ point as shown in Fig. 7(b) by the accidental degeneracy caused by the particular choice of device parameters and then repel each other as in (c) and (d). The two dispersion curves are linear in the vicitiny of the Γ point in Fig. 7(b) as expected from Eq. (54).

Here, let us recall that the singular behavior in Fig. 7(b) was caused by a particular combination of symmetries of the two resonance states. In the present case, they are *A*_{1} and *B*_{1}. The same behavior is expected for a combination of *A*_{2} and *B*_{2} states, since they have the same odd parity with respect to the *y* coordinate. The remaining two combinations, (*A*_{1}, *B*_{2}) and (*A*_{2}, *B*_{1}), have different *σ _{y}* parity so that the corresponding bands do not mix and they merely cross each other.

Finally, the mode frequency at the Γ point depends on the band width of dispersion curves, which is described by *L*_{1} as Eqs. (12) and (14) show. When the lattice constant is increased, the overlapping of wave functions of the nearest neighbors decreases so that the band width also decreases. To confirm this, we calculated the mode frequencies on the Γ point for various lattice constants. The results are plotted in Fig. 8, which clearly shows that (1) the frequency of the *A*_{1} (*B*_{1}) mode decreases (increases) with increasing lattice constant, which is consistent with the decreasing band width, since the *A*_{1} (*B*_{1}) band is concave-down (up) on the Γ point, and (2) the mode frequencies approach the resonance frequencies of a single unit structure (82.8 GHz and 96.6 GHz), when the lattice constant becomes sufficiently large so that the influence of neighboring unit structures is vanishing. When the lattice constant is approximately equal to 0.576 mm, the two curves in Fig. 8 cross each other, which implies that the two mode frequencies coincide as shown in Fig. 7(b). When the lattice constant is smaller than this value and the order of mode frequencies is reversed, dispersion curves should look like Fig. 7(c) and (d), since the parities of the two modes with respect to the *y* coordinate are the same so that the two modes satisfy all conditions to realize the situation of Fig. 7.

## 4. Conclusion

We reformulated the problem of a steerable left-handed antenna from the view point of structural electromagnetic resonance of the metallic unit structure. We derived simple expressions of dispersion curves based on the tight-binding picture by assuming two resonant states of the unit structure with different parities, and described the steering property qualitatively. Through FDTD calculations, we showed the presence of such resonances, which resulted in the formation of two electromagnetic bands with opposite signs of curvature at the Γ point. When the symmetry of the resonance states satisfies certain conditions, the two dispersion curves are linear and have finite slopes in the vicinity of the Γ point in the case of accidental degeneracy. We also showed that the device performance could be designed by changing the lattice constant of the arrayed unit structures, since the dispersion curves were modified due to the change in the electromagnetic transfer integral.

## Acknowledgments

This study was supported by a Grant-in-Aid for Scientific Research on Innovative Areas from the Japanese Ministry of Education, Culture, Sports and Technology (Grant number 22109007).

## References and links

**1. **V. G. Veselago, “Electrodynamics of substances with simultaneously negative values of sigma and mu,” Sov. Phys. Usp. **10**, 509–514 (1968). [CrossRef]

**2. **J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science **312**, 1780–1782 (2006). [CrossRef] [PubMed]

**3. **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**, 4184–4187 (2000). [CrossRef] [PubMed]

**4. **R. A. Shelby, D. R. Smith, and S. Schultz, “Experimental verification of a negative index of refraction,” Science **292**, 77–79 (2001). [CrossRef] [PubMed]

**5. **D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, “Metamaterial electromagnetic cloak at microwave frequencies,” Science **314**, 977–980 (2006). [CrossRef] [PubMed]

**6. **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**, 195104 (2002). [CrossRef]

**7. **S. A. Ramakrishna and T. M. Grzegorczyk, *Physics and Applications of Negative Refractive Index Materials* (SPIE Press, 2008). [CrossRef]

**8. **S. Matsuzawa, K. Sato, Y. Inoue, and T. Nomura, “W-band steerable composite right/left-handed leaky wave antenna for automotive applications,” IEICE Trans. Electron. **E89-C**, 1337–1344 (2006). [CrossRef]

**9. **A. Grbic and G. V. Eleftheriades, “Experimental verification of backward-wave radiation from a negative refractive index metamaterial,” J. Appl. Phys. **92**, 5930–5935 (2002). [CrossRef]

**10. **C. Caloz and T. Ito, “Application of the transmission line theory of left-handed (LH) materials to the realization of a microstrip LH line,” IEEE-AP-S Int. Symp. Dig. **2**, 412–415 (2002).

**11. **S. Tokoro, K. Kuroda, A. Kawakubo, K. Fujita, and H. Fujinami, “Electronically scanned millimeter-wave radar for pre-crush safety and adaptive cruise control system,” *Proc. IEEE Intelligent Vehicles Symp.*, 304–309 (2003).

**12. **J. D. Joannopoulos, R. D. Meade, and J. N. Winn. *Photonic Crystals: Molding the Flow of Light* (Princeton University Press, Princeton, 1995).

**13. **K. Sakoda, *Optical Properties of Photonic Crystals*, 2nd Ed. (Springer-Verlag, Berlin, 2004).

**14. **T. Ito and K. Sakoda, “Photonic bands of metallic systems. II. Features of surface plasmon polaritons,” Phys. Rev. B **64**, 045117 (2001). [CrossRef]

**15. **A. Taflove, *Computational Electrodynamics* (Artech House, Boston, 1995).

**16. **D. M. Sullivan, *Electromagnetic Simulation Using the FDTD Method* (IEEE Press, Piscataway, 2000). [CrossRef]

**17. **N. Peyghambarian, S. W. Koch, and A. Mysyrowicz, *Introduction to Semiconductor Optics* (Prentice Hall, Englewood Cliffs, 1993) Sec. 2.5.

**18. **T. Inui, Y. Tanabe, and Y. Onodera, *Group Theory and Its Applications in Physics* (Springer, Berlin, 1990).

**19. **P. Yeh, “Electromagnetic propagation in birefringent layered media,” J. Opt. Soc. Am. **69**, 742–756 (1979). [CrossRef]

**20. **C. Caloz, A. Lai, and T. Itoh, “The challenge of homogenization in metamaterials,” N. J. Phys. **7**, 167 (2005). [CrossRef]

**21. **A. Lai, T. Itoh, and C. Caloz, “Composite right/left-handed transmission line metamaterials,” IEEE Microwave Magazine, September issue, 34–50 (2004). [CrossRef]