1. Introduction

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 [1].

Graded matematerials (GMs) [2] 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 [6], field rotator [7] and bend structure [8] 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 [16]. The basis of the devices was a dielectric slab structured with holes of variable relative distances [13] 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 [24] and for realization of focusing lens [25] and graded-index waveguide [26].

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.

 

Fig. 1. (a) Refractive index profile with 2D gradient and (b) its discrete approximation. (c) Implementation of the profile with the GPC.

Download Full Size | PDF

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 x_i and y_j and the cell refractive index, n_ij, is equal to n(x_i,y_j). The ray path through the GRIN medium is governed by equations of Hamiltonian optics [27]. 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 [28]. Here equifrequency contours (EFCs) can be approximated with circles and the GPC can be modeled with isotropic and homogeneous medium locally [29], 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, r_ij, of the ij-cell:

 

Fig. 2. Refractive index calculated from dispersion curves of a PC as a function of normalized frequency Ω = a/λ for the PC with (a) SiO₂ and (b) Si rods. On the left side, (x1), the values for TM mode are shown and on the right side, (x2), the values for TE mode. Straight red lines represent values of refractive index calculated with EMT. Ratios r/a are changed from 0.1 to 0.5: □(r/a = 0.1), ○(r/a = 0.2), ♢(r/a = 0.3), △(r/a = 0.4), ∇(r/a = 0.5).

Download Full Size | PDF

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 SiO₂ (ε = 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 [29] 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 [32]. 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 [33] whereas here we use the MG theory [30, 34, 35] so the effective permittivity ε_plane of uniaxial PC is given as [36]

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 ² π/a ² and L_plane = 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 $n_{\mathrm{EMT}}^{\mathrm{TE}}=\sqrt{{\epsilon }_{\mathrm{plane}}}$ and $n_{\mathrm{EMT}}^{\mathrm{TM}}=\sqrt{{\epsilon }_z}$ . 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. [34]. 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 [36]. 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 = n_ij. 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, f_ij = r ² _ij π/a ². Finally, in the case of TE polarization, the radius r_ij 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. [37]. 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.

 

Fig. 3. FEM simulation results of TE mode for (a) the original lens and (b) the GPC based Luneburg lens, Ω₀ = 0.14 whereas Ω ∈ (0.14,0.2). On the left side, (x1), the z-component of magnetic field is shown and on the right side, (x2), the magnetic field intensity.

Download Full Size | PDF

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 n_min, to maximal value n_max, common characteristic for all cells in GPCs which implement devices is free-space wavelength λ ₀, while wavelength λ = λ ₀/n and the normalized frequencies Ω = a/λ vary from cell to cell. Therefore, the excitation frequency in simulations will be expressed in terms of Ω₀ = a/λ ₀, while Ω will lie between n_minΩ₀ and n_maxΩ₀.

3.1. 2D Luneburg lens

The Luneburg lens [40] 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. [41] whereas the implementation of sliced Luneburg lens with drilled holes was presented within the framework of metamaterials in Ref. [42]. Recently, 2D Luneburg lens has been experimentally realized at microwave frequencies using metamaterials based on non-resonant metallic inclusions [43].

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 (n_min = 1, n_max = 1.41), SiO₂ 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.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 Ω₀, Fig. 4, while the reflection from the original lens is given for comparison. As can be seen, for Ω₀ = 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.

 

Fig. 4. (Reflection (S ₁₁ parameter) from the original (blue) and GPC (red) Luneburg lens.

Download Full Size | PDF

Simulation results for the GPC lens at lower frequencies, Ω₀ = 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.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.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 Ω₀ 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 = n_max, implying the lowest value of λ so this is the place where Ω reaches the Bragg reflection condition firstly. With further increasing of Ω₀, 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.3 is the maximal excitation frequency and Ω ≈ 0.42 is the upper frequency limit for proper work of GPC lens.

 

Fig. 5. FEM simulation results of the magnetic field intensity for TE mode in the GPC based Luneburg lens: (a) Ω₀ = 0.08, (b) Ω₀ = 0.3, (c) Ω₀ = 0.33 and (d) Ω₀ = 0.38.

Download Full Size | PDF

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 Ω₀ ≲ Ω_max/n_max. 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%.

 

Fig. 6. The mapping $z=\exp \left(\frac{w}{C}\right)$ , C = 1, for a design of an electromagnetic beam bend: (a) the square domain in Cartesian coordinates of complex plane w = u + iv and (b) its image, the annular segment in curvilinear coordinates of complex plane z = x + iy. The inner and outer bend radii are exp(u ₁) and exp(u ₂), respectively, whereas π/2 is the bend angle. Red arrows show propagation directions of an electromagnetic wave.

Download Full Size | PDF

3.2. Electromagnetic beam bend

The bend is a device intended for steering and guiding of electromagnetic beams. In Ref. [44] we investigated the bend design and implementation in the framework of transformation optics whereas here we designed it by the conformal mapping [45]. The underlying mapping is

which is also used for a design of transformation optics based devices [46]. 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 n_w = 1, refractive index profile within the bend is [45]

Exhaustive analysis of a Gaussian beam propagation through 2D PCs with slowly varying nonuniformities was given in Ref. [18] while the experimental realization of the bend with structured dielectric slab has been recently reported [16]. Here we investigated validity of the proposed method for TM mode. For the chosen bend geometry refractive index varied from n_min = 1.25 to n_max = 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.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.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 [47]. 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].

 

Fig. 7. FEM simulation results of TM mode for (a) the original bend and (b) the GPC based bend, Ω₀ = 0.11. On the left side, (x1), the z-component of electric field is shown and on the right side, (x2), the electric field intensity.

Download Full Size | PDF

The design based on conformal mappings is valid in the limit of geometrical optics [45] so refractive index should vary slowly over the wavelength [27]. 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.068.

As can be seen from Fig. 8, the transmissions for the both bends are at the same level up to Ω₀ = 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.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.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 Ω₀ ≲ Ω_max/n_max.

 

Fig. 8. Transmission (S ₂₁ parameter) through the original (blue) and GPC (red) electromagnetic beam bend. The small oscillations in transmissions are Fabry-Perot resonances within the bends due to impedance mismatching at the bend entrance and exit.

Download Full Size | PDF

 

Fig. 9. FEM simulation results of the z-component of electric field for TM mode in the GPC based bend: (a) Ω₀ = 0.068 and (b) Ω₀ = 0.15.

Download Full Size | PDF

4. Conclusion

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 Ω₀ ≲ Ω_max/n_max. 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.