Engineering of a refractive index profile is a powerful method for controlling electromagnetic fields. In this paper, we investigate possible realization of isotropic gradient refractive index media at optical frequencies using two-dimensional graded photonic crystals. They consist of dielectric rods with spatially varying radii and can be homogenized in broad frequency range within the lowest band. Here they operate in metamaterial regime, that is, the graded photonic crystals are described with spatially varying effective refractive index so they can be regarded as low-loss and broadband graded dielectric metamaterials. Homogenization of graded photonic crystals is done with Maxwell-Garnett effective medium theory. Based on this theory, the analytical formulas are given for calculations of the rods radii which makes the implementation straightforward. The frequency range where homogenization is valid and where graded photonic crystal based devices work properly is discussed in detail. Numerical simulations of the graded photonic crystal based Luneburg lens and electromagnetic beam bend show that the homogenization based on Maxwell-Garnett theory gives very good results for implementation of devices intended to steer and focus electromagnetic fields.
© 2010 Optical Society of America
Metamaterials are artificial electromagnetic structures consisting of unit cells whose size a is much smaller then the wavelength λ so they can be regarded as media with effective parameters. In most of reported realizations, the ratio Ω = a/λ was not negligible but the effective parameters were still well defined. This transitional regime between the effective medium and the photonic crystal regime is denoted as the metamaterial regime .
Graded matematerials (GMs)  consist of spatially varying unit cells and enable realization of various gradient refractive index (GRIN) media for controlling electromagnetic fields. The GRIN media were realized using GMs with metallic inclusions in the cases of electromagnetic cloak [3–5], carpet cloak , field rotator  and bend structure  or using structured metallic surfaces capable for focusing and guiding of surface waves [9–11].
However, in order to fabricate broadband and low-loss GRIN media at optical frequencies, dielectric media should be used instead of metallic inclusions. This has been recently reported for carpet cloak [12–15] and beam steering and focusing devices . The basis of the devices was a dielectric slab structured with holes of variable relative distances  or with rods/holes of variable radii [15,16] and such slab can be modeled as medium with spatially varying refractive index. This way for realizations of GRIN media suggests that they could be implemented by inhomogeneous photonic crystal (PC) with gradually varying unit cells.
Guiding of light in inhomogeneous PCs using not the existence of a photonic band-gap but a well designed spatially dependent dispersion was proposed in Refs. [17, 18]. The concept of graded photonic crystals (GPCs) was introduced in Refs. [19, 20] where it was shown how to bend a path of light by two-dimensional (2D) GPC with one-dimensional (1D) lattice gradient. Similar GPCs were then used for a design of focusing and guiding devices [21–23]. The previous studies on 2D GPCs were focused on realizations of particular functions, but they did not answer how to determine the geometry of unit cells in GPCs which would implement desired GRIN media. Also, the 2D GPCs were graded in one direction giving way for realization of GRIN media with a 1D gradient only. Graded subwavelength gratings could be also regarded as 1D GPCs and they were used for for achieving highly efficient diffractive components  and for realization of focusing lens  and graded-index waveguide .
In this paper we give a general procedure for implementation of specified isotropic GRIN media with 2D GPCs. The GPCs are described as effective media and they operate in metamaterial regime. The effective refractive index of the GPCs is controlled by varying the radii of dielectric rods in both directions which enables realization of GRIN media with a 2D gradient. This extends the scope of applicability of GPCs and enables realizations of more complex GRIN devices. Homogenization of GPCs is approximately done with Maxwell-Garnett (MG) effective medium theory (EMT). Using this approximation, we give analytical formulas to determine the rods radii. The GPCs are low-loss and enable implementation at optical frequencies as well. The proposed method is verified using numerical simulations of GPC based Luneburg lens and electromagnetic beam bend where we focus on investigation of the frequency range where the homogenization procedure is applicable.
2. Graded photonic crystals in metamaterial regime
Consider a 2D device with a GRIN medium whose refractive index profile n(x,y) in xy-plane is shown in Fig. 1(a). In order to realize this profile, it is firstly approximated with discrete one as shown in Fig. 1(b). Each ij-cell is a square of side a, the coordinates of its central point are xi and yj and the cell refractive index, nij, is equal to n(xi,yj). The ray path through the GRIN medium is governed by equations of Hamiltonian optics . Hamiltonian for the ij-cell is a plane wave dispersion
where c is speed of light in vacuum and k is modulus of wave vector in xy-plane.
Implementation of the refractive index profile from Fig. 1(b) requires use of many materials with appropriate discrete values of refractive index. If the refractive index distribution is slowly varying, simpler realization can be achieved by GPC with square lattice whose cross section in the xy-plane is shown in Fig. 1(c). Dielectric rods of only one material are used while their radii are spatially varying. εrod is the permittivity of dielectric rods and εhost is the permittivity of the host medium. The field propagation can be treated by equations of Hamiltonian optics [17,18] if the field can be approximated with a plane wave locally. This is fullfiled in the lowest band except in the vicinity of the band edge . Here equifrequency contours (EFCs) can be approximated with circles and the GPC can be modeled with isotropic and homogeneous medium locally , whose dispersion is
with (nαeff(k))ij being the frequency dependent effective refractive index obtained from dispersion curves of a PC whose unit cell is the ij-cell whereas α stands for transverse magnetic (TM, magnetic field is in xy-plane) and transverse electric mode (TE, electric field is in xy-plane). Equations (1) and (2) have to be equivalent in order that the GPC from Fig. 1(c) would implement GRIN medium from Fig. 1(b), which gives the condition to determine a radius of rod, rij, of the ij-cell:
The ratio between unit cell size a and wavelength λ, Ω = a/λ, (introduced at the beginning of the paper) will be hereafter denoted as normalized frequency. In this paper we studied GPCs with SiO2 (ε = 4.5) and Si (ε = 11.8) rods whose EFCs can be approximated with circles locally up to even Ωmax ≈ 0.44 in the lowest band. This means that the GPCs can be homogenized  and work in metamaterial regime up to Ωmax which is a decisive reason for broadband work of GPC based devices. In Fig. 2 are shown calculated values nαeff(k) (blue lines) of a PC for five different values r/a. The dispersion curves of the PC were calculated using COMSOL Multiphysics by considering single unit cell of the PC with periodic boundary conditions. The eigenvalue problem was then solved numerically while the effective refractive index was calculated from Eq. (2). Due to spatial dispersion, the effective refractive index starts to differ from the value in the long-wavelength limit nαeff. But for the normalized frequencies Ω ≲ 0.25, we can adopt the approximation nαeff(k) ≈ nαeff, so Eq. (3) reads
In the long-wavelength limit, PCs can be regarded as effective homogeneous media [30–32] with a diagonal effective dielectric tensor in principal set of axes. Denote principal dielectric permittivities with εκ(κ = x,y,z). In a general case, the PCs are biaxial, but in the considered case of circular cilynders in a square lattice, they are uniaxial, εx = εy = εplane, that is, they are isotropic in xy-plane . In this paper we use analytical formulas for homogenization of PCs in the long-wavelength limit since it significantly simplifies design of GPC based devices. Homogenization of periodic structures using analytical formula has been done long ago  whereas here we use the MG theory [30, 34, 35] so the effective permittivity εplane of uniaxial PC is given as 
while εz is obtained using the exact formula, that is, as the volume average of εrod and εhost
Here f stands for the filling fraction of rods, f = r 2 π/a 2 and Lplane = 1/2 is the depolarizing factor of rods in the long-wavelength limit for xy-plane. The effective refractive indices based on EMT theory are and . Their values are represented with straight, red lines in Fig. 2. As can be seen, in the long-wavelength limit, there is a complete agreement with values obtained from dispersion curves for TM mode while for TE mode, acceptable agreement is observed for r/a < 0.4. This is in accordance with the results reported in Ref. . For TE mode, the electric field is normal to the rods in PCs and the main assumption in the derivation of MG theory is that a rod is placed in the local field of homogeneously polarized matter . For low filling fraction, this is a good approximation, but with increasing radius of rods, the field within the rod becomes more perturbed and the approximation of homogeneously polarized environment fails. This is the reason for the significant difference between the effective index for TE mode obtained from MG theory and the one calculated from dispersion curves for r/a > 0.4 even in the long-wavelength limit. In the case of TM mode, the electric field is parallel to the rods so the PC is homogenously polarized even for large filling fractions and field averaging give exact value in the long-wavelength limit.
Since the effective index nαeff can be expressed using EMT, the condition from Eq. (4) now reads (nαEMT)ij = nij. Radii of rods in the GPC from Fig. 1 (c) can be determined by putting this condition in Eqs. (5) and (6) with spatially dependent filling fraction, fij = r 2 ij π/a 2. Finally, in the case of TE polarization, the radius rij is obtained from Eq. (5) as
whereas for TM mode, the radius is obtained from Eq. (6) as
In the previous discussion we neglected frequency dispersion of rods permittivity so Eqs. (7) and (8) should be used in frequency range where this approximation is valid. The use of only dielectric inclusions makes possible realization of low-loss devices but limits values of effective refractive index achievable with GPCs. The lower bound is 1 while the upper bounds are determined with a maximum allowed filling fraction. In the case of TM mode, the maximum filling fraction is 0.78 because of the touching rods whereas in the case of TE mode, the maximum filling fraction is 0.5 due to the condition r/a < 0.4.
The analytical formulas [Eq. (7) and (8)] greatly simplify the design compared to that based on extracting effective index from numerical simulations of unit cells and adjusting their geometries to fit the specified refractive index. The main assumption was the approximation nαeff(k) ≈ nαEMT for Ω ≲ 0.25 which is consistent with the result reported in Ref. . The approximation is applicable due to approximately linear dispersion curves of PC in the lowest band. In the next section, it will be shown that this approximation could give good results in implementation of GRIN media even for frequencies Ω ≳ 0.25, although here the difference between nαeff(k) and nαEMT is more pronounced. By calculating nαEMT using extended MG theories [38, 39], it would be possible to completely match frequency dependent refractive index nαeff(k) in Fig. 2, not only its value in the long-wavelength limit but the effective index nαEMT would be then Ω dependent. On the other hand, Eqs. (7) and (8), enable frequency independent design applicable in broad frequency range.
3. Device implementation by graded photonic crystals in metamaterial regime
The proposed method was validated by numerical simulations of GPC based Luneburg lens and electromagnetic beam bend with finite element method (FEM) based software COMSOL Multiphysics. Since refractive index in the devices changes from minimal nmin, to maximal value nmax, common characteristic for all cells in GPCs which implement devices is free-space wavelength λ 0, while wavelength λ = λ 0/n and the normalized frequencies Ω = a/λ vary from cell to cell. Therefore, the excitation frequency in simulations will be expressed in terms of Ω0 = a/λ 0, while Ω will lie between nminΩ0 and nmaxΩ0.
3.1. 2D Luneburg lens
The Luneburg lens  is a spherically symmetric lens which focuses incoming electromagnetic field from any direction to the point at the opposite lens side or transforms a radiation of a point source at the lens surface into parallel rays at the opposite lens side. The lens is very attractive as a receiver and transmitter in many antenna applications and as a field concentrator for focusing. Standard realization is based on spherical dielectric shells with constant permittivity. Using of PCs for optical devices such as lens was proposed in Ref.  whereas the implementation of sliced Luneburg lens with drilled holes was presented within the framework of metamaterials in Ref. . Recently, 2D Luneburg lens has been experimentally realized at microwave frequencies using metamaterials based on non-resonant metallic inclusions .
Here we considered the implementation of 2D cylindrically symmetric Luneburg lens with GPCs as an example for testing the proposed method for TE polarization. The lens refractive index is given by
where R is the lens radius and ρ is modulus of radius vector of the lens internal points. In order to implement this permittivity profile (nmin = 1, nmax = 1.41), SiO2 rods (ε = 4.5) were used whereas the host medium was vacuum. Radii of rods were calculated using Eq. (7). The simulation results for the original lens with refractive index profile given by Eq. (9) are shown in Fig. 3(a) whereas the results for the GPC lens are given in Fig. 3(b). Excitation frequency in the simulated case was Ω0 = 0.14 while Ω ∈ (0.14,0.2). As can be seen the GPC lens works as well as the original one, the plane wave from the left side is focused onto the opposite, right side, while the field within the lens behaves as a plane wave locally which indicates that the effective medium approximation is valid. To give exact criterion for the description of the lens performances and to suggest possible way for experimental testing, we calculated reflection in front of the GPC lens as a function of excitation frequency Ω0, Fig. 4, while the reflection from the original lens is given for comparison. As can be seen, for Ω0 = 0.14, the reflections from both lenses are negligible meaning that the incident field does not resolve periodic structures of the GPC lens perceiving it as a locally homogeneous medium such as the original lens.
Simulation results for the GPC lens at lower frequencies, Ω0 = 0.08, Ω ∈ (0.08,0.11), are shown in Fig. 5(a). As can be seen, the lens works fine while the reflection is low, Fig. 4. However, the focus became wider, it was moved to the lens interior and field intensities in the focus were decreased. These lower performances are because the original lens was designed in the limit of geometrical optics. Nevertheless, simulations showed that the lens worked well even for Ω0 = 0.067, far beyond the geometrical optics when free-space wavelength was equal to the lens radius. It should be emphasized that existence of the lower frequency limit is not consequence of the implementation with GPCs.
The GPC lens performances at higher frequency can be explored through the simulation for Ω0 = 0.3, Ω ∈ (0.3,0.42), where the reflection is still low and at the same level as for the original lens, Fig. 4. As can be seen from Fig. 5(b), the GPC lens works very fine for this frequency, focal point is narrower, field intensities are greater, so focusing is even better since the work at higher frequencies is closer to the limit of geometrical optics. But for Ω ≳ 0.25, the difference between the effective index based on MG theory and the effective index calculated from dispersion curves of PC results in slight changes of refractive index profile within the GPC lens so the focus point is moved toward the lens interior. Therefore, for the frequencies Ω ≳ 0.25, implementation of the Luneburg lens with GPCs is still possible since it is not strongly sensitive to slight changes of refractive index distribution but for the precise control, the effective index should be calculated from dispersion curves of PCs.
Further increasing of excitation frequency Ω0 above 0.3 leads to strong reflection from the GPC lens, Fig. 4. Since the reflection from the original lens is still negligible, the increased reflection means that the incident field does not perceive the GPC lens as a locally homogeneous medium anymore. This can be also observed as a standing wave pattern in a front of and within the central part of the GPC lens, Fig. 5(c), as a result of Bragg reflections from the central cells. The effective refractive index is the largest in the center of the lens, n = nmax, implying the lowest value of λ so this is the place where Ω reaches the Bragg reflection condition firstly. With further increasing of Ω0, the Bragg reflection condition becomes fulfilled in parts of the GPC lens with lower value of effective refractive index, so Bragg reflections start to appear from the side cells in the GPC lens as well, Fig. 5(d). Therefore, based on calculation of reflection from the GPC lens, we conclude that Ω0 ≈ 0.3 is the maximal excitation frequency and Ω ≈ 0.42 is the upper frequency limit for proper work of GPC lens.
The maximal excitation frequency for proper work of the GPC lens can be determined with simple approximative formula. The upper frequency limit Ω ≈ 0.42 after which Bragg reflections appear is very close to Ωmax ≈ 0.44, the maximal normalized frequency for the homogenization and work in metamaterial regime for PCs. This means that Bragg reflections firstly appear when Ω starts to overcome Ωmax within the part of a GPC based device with the highest effective refractive index. For the lens, this is the central part. Therefore, for the proper work of the device, the excitation frequency have to satisfy the approximative condition Ω0 ≲ Ωmax/nmax. The advantage of this formula is because Ωmax can be simply obtained by calculating EFC for unit cell which implements the highest effective refractive index.
It should be noted that the upper frequency limit depends on resolution of a discrete refractive index profile. When the discrete profile approximates the original one sufficiently well, further decreasing of unit cell side will not significantly improve device performances but it will enable work at higher frequencies due to scaling law of PCs. Of course, the price will be more complicated fabrication of smaller inclusions in GPCs. In previous examples, the lens was implemented with 15 unit cells per radius, which implied that the change of refractive index per unit cells was 0.027. Simulations of the lens (not shown here) implemented with 10 (approximately 33% decreased) unit cells per radius and the change of refractive index per unit cells of 0.041 showed that it worked fine, but the upper frequency limit was decreased by approximately 33%.
3.2. Electromagnetic beam bend
The bend is a device intended for steering and guiding of electromagnetic beams. In Ref.  we investigated the bend design and implementation in the framework of transformation optics whereas here we designed it by the conformal mapping . The underlying mapping is
which is also used for a design of transformation optics based devices . By this mapping, square domain from straight Cartesian grid in w-complex plane, w = u + iv, is transformed to annular segment in polar coordinates of z-complex plane, z = x + iy, as shown in Fig. 6. The constant C can be used for adjusting dielectric profile within the bend. If the square domain is empty that is, if refractive index within square domain is nw = 1, refractive index profile within the bend is 
Exhaustive analysis of a Gaussian beam propagation through 2D PCs with slowly varying nonuniformities was given in Ref.  while the experimental realization of the bend with structured dielectric slab has been recently reported . Here we investigated validity of the proposed method for TM mode. For the chosen bend geometry refractive index varied from nmin = 1.25 to nmax = 2.86. This is a significantly larger variation of the refractive index compared to the Luneburg lens so a denser distribution of rods was required in order to obtain good behaviour. The bend was implemented with Si (ε = 11.8) rods in vacuum background, Eq. (8) was used to determine radii of rods while C = 1. The simulation results for the original and GPC based bend are shown in Fig. 7(a) and 7(b), respectively, with excitation frequency Ω0 = 0.11 while Ω∈ (0.137,0.32). As can be seen, the field distribution in the GPC bend matches the field in the original one and the field perceives the GPC bend as an effective medium. In order to exactly compare performances of GPC bend and the original one, we calculated transmissions through both of them which could be also appropriate for experimental testing. The calculation results are shown in Fig. 8 and for Ω0 = 0.11 both bends have the same transmission.
However, both bends suffer from the impedance mismatching at the entrance and the exit which leads to the reflection of the incident field which results in decreased transmission, Fig. 8, while the standing wave appears in front of the bend as well as within it, Figs. 7(a2) and 7(b2). The reflection could be decreased by designing appropriate antireflection coatings or by using numerical conformal transformations for a design which results in reflectionless devices for TM mode . Implementation of the later approach requires permittivity below one within some parts of the devices so metallo-dielectrics structures have to be used. Also, the question is if it is always possible to find appropriate conformal/quasiconformal transformation for a device design. In cases where it is not, ray-tracing method based on Eikonal equation could be used instead [8, 16, 48].
The design based on conformal mappings is valid in the limit of geometrical optics  so refractive index should vary slowly over the wavelength . Although the gradient of refractive index in the bend was 3.5 times larger than in the Luneburg lens, the simulations showed that the bending effect worked even in the cases where wavelength was comparable to the bend dimensions while the transmissions stay at the same level for both bends, Fig. 8. However, a beam spreading at the lower frequencies was more distinct so the beam exceeded the bend width significantly causing distortion of wave front. Such case is shown in Fig. 9(a) for Ω0 = 0.068.
As can be seen from Fig. 8, the transmissions for the both bends are at the same level up to Ω0 = 0.15 (Ω is between 0.187 and 0.43). The simulation results for this excitation frequency are shown in Fig. 9(b). The bending effect works well, but there is a little phase delay (about half a wavelength) in the GPC bend compared to the original one (simulation not shown here). This indicates a slight difference between original distribution of refractive index and the one implemented with GPC as a result of applied homogenization based on EMT. The bending effect showed robustness to this difference but in cases where refractive index have to be precisely controlled, it should be calculated from PC dispersion curves.
Beyond the frequency Ω0 ≈ 0.15 transmission for the GPC bend decreases abruptly while it stays at the same level for the original bend, Fig. 8. This means that the incoming field does not perceive the GPC bend as an effective medium anymore. Bragg reflections start to appear from the cells along the inner bend edge because the refractive index is the highest and the wavelength is the lowest here and this is the place where Bragg condition is satisfied firstly. Since the bending effect does not work anymore, Ω0 ≈ 0.15 is maximal allowed excitation frequency while Ω ≈ 0.43 is the upper frequency limit for proper work of GPC bend. As in the case of the Luneburg lens, the upper frequency limit Ω ≈ 0.43 is very close to Ωmax. Therefore, achieving of that limit can be again interpreted as a result of increasing of Ω above Ωmax for the cells with the highest effective refractive index in GPC. As a result, the maximal allowed excitation frequency can be again determined using the approximative formula Ω0 ≲ Ωmax/nmax.
In a summary, the proposed method for implementation of isotropic GRIN media using GPCs in metamaterial regime was validated through numerical simulations of the Luneburg lens and the beam bend. It was shown that they can be implemented with dielectric rods of only one material by changing their radii. This makes possible fabrication of low-loss devices at optical frequencies. Simplicity of the proposed method is a result of approximately linear dispersion of PCs for wide frequency range in the lowest band. This enables calculation of effective parameters in the framework of MG theory and analytical formulas for determination of rods radii. As a conclusion, the GPC studied in this paper can be regarded as a class of graded dielectric metamaterials.
The proposed method enables work in broad frequency range in the lowest band. For normalized frequencies Ω ≲ 0.25, the GPC based devices worked as well as the original GRIN devices so here the proposed method for implementation gave excellent results. At higher frequencies, Ω ≳ 0.25, homogenization of PCs is still possible since EFCs can be approximated with circles. Although effective parameters of PCs then start to differ from the values obtained from MG theory, in the examples of the Luneburg lens and the bend structure it was shown that the method gave good results and it was very robust. The conclusion is that the method is applicable for implementation of devices intended to guide, steer and focus electromagnetic fields, which are insensitive to a slight variation of refractive index. But, for Ω ≳ 0.25, in the application where it is necessary to control refractive index very precisely, it have to be calculated from dispersion curves.
The upper frequency limit for proper work of GPC based devices was detected with appearance of Bragg reflections from GPC devices. This happens when normalized frequency Ω in the parts of the devices with the highest value of refractive index starts to overcome Ωmax, the maximal normalized frequency for the homogenization and work in metamaterial regime for GPCs. As a result, the maximal excitation frequency for proper work of the GPC devices can be determined using approximative formula Ω0 ≲ Ωmax/nmax. Although the original devices with GRIN medium were designed in the limit of geometrical optics, simulations showed that the GPC based bend and lens worked very well outside the limit of geometrical optics.
This work is supported by the Serbian Ministry of Science under Project No. 141047. G.I. acknowledges support from ORSAS in the U. K. and the University of Leeds. R.G. acknowledges support from EU FP7 projects Nanocharm and NIMNIL. K.H. is grateful to the Austrian NILmeta-NILAustria project from FFG for partial support. We thank Johann Messner from the Linz Supercomputer Center for technical support.
References and links
1. D. R. Smith, D. C. Vier, T. Koschny, and C. M. Soukoulis, “Electromagnetic parameter retrieval from inhomogeneous metamaterials,” Phys. Rev. E 71, 036617 (2005). [CrossRef]
2. D. R. Smith, J. J. Mock, A. F. Starr, and D. Schurig, “Gradient index metamaterials,” Phys. Rev. E 71, 036609 (2005). [CrossRef]
4. 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]
5. W. Cai, U. K. Chettiar, A. V. Kildishev, and V. M. Shalaev, “Optical cloaking with metamaterials,” Nat. Photonics 1, 224–227 (2007). [CrossRef]
7. H. Chen, B. Hou, S. Chen, X. Ao, W. Wen, and C. T. Chan, “Design and experimental realization of a broadband transformation media field rotator at microwave frequencies,” Phys. Rev. Lett. 102, 183903 (2009). [CrossRef]
10. A. O. Pinchuk and G. C. Schatz, “Metamaterials with gradient negative index of refraction,” J. Opt. Soc. Am. A 24, A39–A44 (2007). [CrossRef]
11. B. K. Juluri, S. Chin, S. Lin, T. R. Walker, L. Jensen, and T. J. Huang, “Propagation of designer surface plasmons in structured conductor surfaces with parabolic gradient index,” Opt. Express 17, 2997–3006 (2009). [CrossRef]
13. J. Valentine, J. Li, T. Zentgraf, G. Bartal, and X. Zhang, “An optical cloak made of dielectrics,” Nat. Mater. 8, 569–571 (2009). [CrossRef]
14. L. H. Gabrielli, J. Cardenas, C. B. Poitras, and M. Lipson, “Silicon nanostructures cloak operating at optical frequencies,” Nat. Photonics 3, 461–463 (2009). [CrossRef]
15. J. H. Lee, J. Blair, V. A. Tamma, Q. Wu, S. J. Rhee, C. J. Summers, and W. Park, “Direct visualization of optical frequency invisibility cloak based on silicon nanorod array,” Opt. Express 17, 12922–12928 (2009). [CrossRef]
16. Z. L. Mei, J. Bai, and T. J. Cui, “Gradient index metamaterials realized by drilling hole arrays,” J. Phys. D Appl. Phys. 43, 055404 (2010). [CrossRef]
17. P. S. J. Russell and T. A. Birks, “Hamiltonian optics of nonuniform photonic crystals,” J. Lightwave Technol. 17, 1982 (1999). [CrossRef]
18. Y. Jiao, S. Fan, and D. A. B. Miller, “Designing for beam propagation in periodic and nonperiodic photonic nanostructures: Extended Hamiltonian method,” Phys. Rev. E 70, 036612 (2004). [CrossRef]
20. E. Centeno, D. Cassagne, and J.-P. Albert, “Mirage and superbending effect in two-dimensional graded photonic crystals,” Phys. Rev. B 73, 235119 (2006). [CrossRef]
22. F. S. Roux and I. De Leon, “Planar photonic crystal gradient index lens, simulated with a finite difference time domain method,” Phys. Rev. B 74, 113103 (2006). [CrossRef]
24. S. Astilean, P. Lalanne, P. Chavel, E. Cambril, and H. Launois, “High-efficiency subwavelength diffractive element patterned in a high-refractive-index material for 633 nm,” Opt. Lett. 23, 552–554 (1998). [CrossRef]
25. U. Levy, M. Abashin, K. Ikeda, A. Krishnamoorthy, J. Cunningham, and Y. Fainman, “Inhomogeneous dielectric metamaterials with space-variant polarizability,” Phys. Rev. Lett. 98, 243901 (2007). [CrossRef]
26. U. Levy, M. Nezhad, H.-C. Kim, C.-H. Tsai, L. Pang, and Y. Fainman, “Implementation of a graded-index medium by use of subwavelength structures with graded fill factor,” J. Opt. Soc. Am. A 22, 724–733 (2005). [CrossRef]
27. Y. A. Kravtsov and Y. I. Orlov, Geometrical Optics of Inhomogeneous Media, (Springer-Verlag, 1990).
28. W. Smigaj and B. Gralak, “Validity of the effective-medium approximation of photonic crystals,” Phys. Rev. B 77, 235445 (2008). [CrossRef]
29. P. A. Belov and C. R. Simovski, “Homogenization of electromagnetic crystals formed by uniaxial resonant scatterers,” Phys. Rev. E 72, 026615 (2005). [CrossRef]
30. S. Datta, C. T. Chan, K. M. Ho, and C. M. Soukoulis, “Effective dielectric constant of periodic composite structures,” Phys. Rev. B 48, 14936–14943 (1993). [CrossRef]
32. P. Halevi, A. A. Krokhin, and J. Arriaga, “Photonic Crystal Optics and Homogenization of 2D Periodic Composites,” Phys. Rev. Lett. 82, 719–722 (1999). [CrossRef]
33. L. Lewin, “The electrical constants of a material loaded with spherical particles,” Proc. Inst. Elec. Eng. 94, 65–68 (1947).
34. M. J. A. De Dood, E. Snoeks, A. Moroz, and A. Polman, “Design and optimization of 2D photonic crystal waveguides based on silicon,” Opt. Quantum Electron. 34, 145–159.
35. A. Kirchner, K. Busch, and C. M. Soukoulis, “Transport properties of random arrays of dielectric cylinders,” Phys. Rev. B 57, 277–288 (1998). [CrossRef]
36. A. Sihvola, Electromagnetic mixing formulas and applications, (The Institution of Electrical Engineers, London, United Kingdom, 1999). [CrossRef]
37. W. G. Egan and D. E. Aspnes, “Finite-wavelength effects in composite media,” Phys. Rev. B 26, 5313–5320 (1982). [CrossRef]
38. W. T. Doyle, “Optical properties of a suspension of metal spheres,” Phys. Rev. B 39, 9852–9858 (1989). [CrossRef]
39. R. Ruppin, “Evaluation of extended Maxwell-Garnett theories,” Opt. Commun. 182, 273–279 (2000). [CrossRef]
40. R. K. Lüneburg, The mathematical theory of optics, (University of California Press, Los Angeles, CA, 1944).
41. P. Halevi, A. A. Krokhin, and J. Arriaga, “Photonic crystals as optical components,” Appl. Phys. Lett. 75, 2725–2727 (1999). [CrossRef]
42. G. Zouganelis and D. Budimir, “Effective dielectric constant and design of sliced Luneberg lens,” Microwave Opt. Technol. Let. 49, 2332–2337 (2007). [CrossRef]
43. Q. Cheng, H. F. Ma, and T. J. Cui, “Broadband planar Luneburg lens based on complementary metamaterials,” Appl. Phys. Lett. 95, 181901 (2009). [CrossRef]
44. B. Vasić, G. Isić, R. Gajić, and K. Hingerl, “Coordinate transformation based design of confined metamaterial structures,” Phys. Rev. B 79, 085103 (2009). [CrossRef]
46. S. Han, Y. Xiong, D. Genov, Z. Liu, G. Bartal, and X. Zhang, “Ray optics at a deep-subwavelength scale: a transformation optics approach,” Nano Lett. 8, 4243–4247 (2008). [CrossRef]
48. Z. L. Mei and T. J. Cui, “Arbitrary bending of electromagnetic waves using isotropic materials,” J. Appl. Phys. 105, 104913 (2009). [CrossRef]