## Abstract

We explore two-dimensional triangular lattice photonic crystals composed of air holes in a dielectric background which are subject to a graded-index distribution along the direction transverse to the propagation. The proper choice of the parameters such as the input beam width, gradient coefficient, and the operating frequency allow the realizations of the focusing (lens) and guiding (waveguide) effects upon which more complex optical devices such as couplers can be designed. Numerical results obtained by the finite-difference time-domain and planewave expansion methods validate the application of Gaussian optics within a range of parameters where close agreement between them are observed.

©2007 Optical Society of America

## 1. Introduction

Periodic dielectric materials known as photonic crystals (PCs) have become indispensable design components in many areas of technology [1–4]. Building novel and realistic functionalities based upon PCs opens applications difficult to realize by other means. With this motivation in mind, in the following we explore the concept of graded-index (GRIN) PCs. It is known that in smoothly graded-index media, varying the refractive index around the optical axis can focus or defocus the input beam. The profile of this variation leads to particularly interesting properties when that profile is quadratic. Slowly varying the refractive index in the direction transverse to the propagation has been exploited for coupling and mode matching in fiber optics [5–7]. In PCs, the band-gap and the super-collimation effects are the two main mechanisms employed to guide light. The first relies on perturbing the specific cells in a PC, within a line or within a linear chain of tightly coupled cells. As a result, a PC waveguide (PCW) or coupled-resonator optical-waveguide is obtained [8–11]; however, one may not need a geometric perturbation for super-collimation, if the operating frequency is outside of the band gap [12]. The equifrequency (EF) contours can be used to infer information about the direction of light propagation inside the PC. The light propagation direction is perpendicular to these contours due to the *ν _{g}* = ∇

*(*

_{k}ω*k*), with

*ν*the group velocity. The contour shapes can be simple or complex depending on the type of PC, polarization, and frequency. If the contour is flat, then self-collimation may take place [12]. This property has been reported earlier and device applications utilizing self-collimation have been proposed [13, 14]. In our case, self-collimation is obtained at two different frequency regimes, below and above the bandgap. Mainly, focusing and collimation effects are governed by graded index medium (details are shown later) which is a new approach.

_{g}Recently, Centeno *et al*. studied graded photonic crystals [15, 16]. In their work index variation (modulation) is obtained by linearly perturbing the lattice spacing (increasing) in the *y* direction. One main difference between our work and theirs in addition to the way of obtaining index variation (parabolic in our case) is the direction of the wave propagation and index variation direction. They propagate the input Gaussian beam along the index variation direction. They have half of the graded index medium obtained by linear index variation. The index variation direction is transverse to the Gaussian beam propagation direction in our case. Also, the index profile is symmetric with respect to the center of the PC which is not the case in previously published results. There are other efforts to engineer PC for focusing and dispersion engineering goals [17–20].

Here, we propose a novel way to achieve focusing and self-guiding of light using two-dimensional (2D) PCs with perforated holes in a high-index dielectric background by merging the idea of graded-index waveguides and PCWs. The demand for compact focusing and self-collimated guiding is high. Usually lenses used for focusing are bulky and susceptible to misalignment. For the self-collimated behavior, one may need to depend on the nonlinear mechanism like Kerr nonlinearity to beat the diffraction. However, the idea presented in this study to achieve the aforementioned light behaviors is free from the nonlinear effects and can be implemented in an integrated manner. The index variation within the medium occupying the holes is taken to be quadratic with respect to the central section of the PC, i.e., *y*-direction (ΓM) and light propagates in *x*-direction (ΓK) as shown in Fig. 1. Numerical analysis is carried out by the 2D finite-difference time-domain method (FDTD) and the planewave expansion technique (PWM), and they are verified within a range of validity of Gaussian optics [21, 22]. It is important to note that as long as the index variation is small within one wavelength, Gaussian optics is valid. Gaussian optics analysis is employed for two main reasons. One is for convenience; there is a wealth of analytical expressions enabling a full understanding of the beam behavior. The first reason would be merely a curiosity were it not for the second, namely, that Gaussian optics predicts stable, well-behaved beams in quadratically graded media. That said, we show that similarly stable, well-behaved beams exist in GRIN PCs. The effective index method is used to fit the results of GRIN PC (from FDTD) to standard Gaussian optics. The GRIN PC, however, enables frequency regimes where Gaussian optics may often not in practice be applicable. The advantages of photonic crystals further include possibly more flexible geometrical manipulations (i.e., guiding in various directions, coupling to other devices, etc.) by standard fabrication methods. In addition, the concept will enable integration of various PC structures on a single substrate and may further facilitate the coupling or integration of various types of guiding structures such as conventional integrated optics, photonic-crystal waveguides, freespace optics, etc. This work is an initialization step to integrate graded index photonic crystals in device applications. For example, one can use focusing properties of GRIN PC to couple spatially wide input beams to photonic crystal waveguides (mode matching). This aspect is out of the scope of the study and target of the future investigation. One should again note that PCs provide GRIN medium profiles by geometrical manipulations which could be difficult to achieve with any other methods, e.g., ion implantation in fiber optics.

Other effective index variations are clearly possible; however, in many cases the leading contribution to the index variation is naturally quadratic (e.g., by ion implantation) [23]. But a leading motivation is that stable Gaussian beams exist in quadratic media, and at least to lowest order, a Gaussian variation may be a good approximation to structures containing high-order variations as well [24].

## 2. Design and analysis of graded index photonic crystals

There are many options to implement the graded-index variation within the PC, for example by changing the radii of air holes. Another way is to modulate the refractive index of the holes. The realization of such index variation on top of the crystalline structure of PC as shown in Fig. 1(a) can be obtained either by varying radius (r) in each row [Fig. 1(b)] or modulating hole refractive indices [Fig. 1(c)]. A schematic of the proposed structure is shown in Fig. 2 along with the desired index variation across a *y* section. The index profile within the holes is parabolic and has its maximum value at the center of the PC. It is assumed to be of the form of *n*
^{2}(*y*) = *n*
^{2}
_{0}⌊1-*α*
^{2}(*y*-*y*
_{0})^{2}⌋, where *α* is the gradient coefficient and *n*
_{0} is the refractive index at *y* = *y*
_{0}. For a case study *n*
_{0} and α values are taken as 3 and 0.09*a*
^{-1}, respectively. The underlying PC has dielectric background with *ε* = 12, and hole-diameter to the lattice-constant ratio is taken as 0.4*a*. The polarization is assumed to be TE (magnetic field is along the air holes), i.e., the electric field is in the periodic plane. Each unit cell is divided into 30*×*30 discretization grid cells.

Gaussian beam propagation through lenslike media has been studied extensively [24, 25]. The analytical framework to study beam propagation through such a medium is Gaussian optics. The details of the following formulation are provided for the completeness of the study. There, one defines the complex beam parameter,

where *R*(*x*) is the wave-front curvature, *ω*(*x*) is the beam spot-size, and *λ*
_{0} is the free-space wavelength [25]. The output beam parameter at *x* = *d* is related to the input beam parameter by the ABCD law via the following relation,

where *r*
_{2} is the ray position and *r*
_{2}́ is the ray slope in the output plane and the corresponding parameters in the input plane are represented by *r*
_{1} and *r*
_{1}́, respectively. The matrix is

From the above matrix, the oscillatory behavior of the output is obvious with periodicity of 2*π*/*α*. The transformation law as stated in Eq. (2) for a Gaussian beam implies

Based on an analysis of Eqs. (3) and (4), there are three distinct possibilities depending on the width of the input pulse: the field may diverge, converge, or maintain a constant beam width (stays collimated) as we keep *α* constant. The oscillatory behavior for the first two cases exhibits a period that is inversely proportional to *α*.

A special circumstance occurs when the initial beam is a planewave, i.e., *R*(*x* = 0) → ∞. Then, the input and output beam parameters are

The condition to have a collimated beam is that *R*
_{out} → ∞ regardless of the *x*-coordinate (propagation direction). In that case,

Then, after some algebra, the output beam width is

It can be seen that frequency and beam width are inversely proportional to *w*
_{0}.

To study the second special case, the imaginary part of Eq. (6) is equated with the complex beam parameter of the output field with the assumption that the phase fronts are planar at the foci,

where *n _{x}* is the index of the medium which is equal to

*n*

_{0}for the conventional approaches but it is different with PCs (the value of the

*n*is approximated to replace the PC with an effective index). Then,

_{x}where ${q}_{\mathrm{in}}=i\phantom{\rule{.2em}{0ex}}\frac{\pi {n}_{0}{w}_{0}^{2}}{{\lambda}_{0}}=i{q}_{i}^{\prime}$. Inserting this value into the above equation gives

From this equation, we conclude that the variation of the width of the propagating field through the graded-index medium is controlled by the parameters *w*
_{0}, *α*, *n*
_{0}, *n _{x}* and is frequency dependent.

The general form of the field in GRIN media is the solution of the paraxial wave equation, and it can be written as

$$\phantom{\rule{.1em}{0ex}}={H}_{0}\frac{{w}_{0}}{w\left(x\right)}\mathrm{exp}\left(-i\phantom{\rule{.2em}{0ex}}\frac{\pi {n}_{0}{y}^{2}}{{\lambda}_{0}R\left(x\right)}\right)\mathrm{exp}\left(-\frac{{y}^{2}}{{w}^{2}\left(x\right)}\right)\phantom{\rule{8.6em}{0ex}}$$

where *H*
_{0} is a constant. It describes the amplitude variation of the magnetic field. The field amplitude and its width are inversely proportional (power must be constant and is independent of *x*). The above analysis will provide insight as we move to study the GRIN PC.

As a first step, the GRIN PC shown in Fig. 2 is studied under different input beam widths *w*
_{0}. The source (modulated Gaussian pulse) is located in front of the GRIN PC centered at the axis of the symmetry (*y* = *y*
_{0}). We mean by the pulse throughout the manuscript is the modulated Gaussian at some center frequency with Δ*f*/*f*
_{0} ≅ 0.067, where *f*
_{0} is the center frequency and Δ*f* is the full-width at half-maximum value. The center frequency is varied to trace the frequency axis in the dispersion diagram. But the bandwidth of the modulated Gaussian pulse is kept constant, only the spatial width (FWHM) along *y*-direction is changed to study the beam width effect on the light propagation through GRIN PC. Three distinct behaviors are observed. For small and large values of *w*
_{0} the pulse diverges and converges, respectively. However, for a specific intermediate value of *w*
_{0}, the field exhibits constant output pulse width. Figures 3(a) and 3(b) show the steady-state field variations of the self-collimated and focusing behavior of GRIN PC at frequencies *a*/*λ* = 0.12 and *a*/*λ* = 0.24, respectively. To investigate the responsible mechanism behind the self-collimation in Fig. 3(a), we plot EF contours in Fig. 4 for the lowest three bands at two different cases (*ε _{a}* = 1.0 and

*ε*= 4.0, where

_{a}*ε*is the hole dielectric index value). It can be seen that the modification of holes refractive index induces local modification of EF curves while keeping their shapes the same. By looking at the EF curve which is circular (dotted circle) in Fig. 4 at

_{a}*a*/

*λ*= 0.12 case, we can attribute the self collimation effect not due to the flat dispersion curve but due to the graded index modulation. Figure 4 will be referred again later to explain other situations observed with GRIN PC.

The beam profile at different propagation distances are plotted in Fig. 5. It can be seen that beam does not spread for the full-width at half-maximum (FWHM) of 2.90 *μ*m [Fig. 5(a)] if *λ* is taken as 1.55 *μ*m. The beam profile for the focusing case is shown in Fig. 5(b) at the input (solid line) and the focal point (dashed line). The FWHM values are 4.22 *μ*m and 0.97 *μ*m, respectively. The spot-size converter ratio defined as *w _{x}* :

*w*

_{0}is 4.35:1. The profiles of the focused beam as well as of the collimated beam closely resemble a Gaussian beam. To explore the characteristics of the GRIN PC and implement it in the device applications (lens and waveguide), we further explored the cases in Figs. 3(a) and 3(b) in Sections 3 and 4.

## 3. Lens design

The converging behavior of the GRIN PC under appropriate conditions enables one to design a lens that focuses the input beam to a focal point [Figs. 3(b) and 5(b)]. There are two key design parameters: *w*
_{0} and *α*. For a range of *w*
_{0} values the GRIN PC brings the input field to focus and the focal length is inversely proportional to *α*. The focal plane predicted by paraxial approximation obeys the formula

Larger *α* value brings the field to a focus in a shorted distance with smaller spot-size.

The beam waist at the focal point *w _{x}* can be adjusted by varying

*w*

_{0}and

*α*. The flexibility to modify the beam waist and adjust the focal point is obviously a great advantage for coupling light between different components of integrated optical systems. A large spot- size ratio can be obtained, which will be of great utility for input and output coupling devices.

## 4. Waveguide design

The other special circumstance happens if the beam width neither diverges nor converges as shown in Fig. 3(a). The condition for this situation can be obtained from Eq. (8) by assuming that the radius of the curvature is very large. At a specific beam size, the diffraction of the Gaussian beam and focusing power of the GRIN PC balance each other. As a result, waveguides with constant beam width (no divergence) can be obtained. Figure 5(a) shows the beam profile at three different positions along the propagation axis. It can be seen that the Gaussian profile of the beam remains unchanged as it travels through the GRIN PC.

In the above two examples, the operating frequency of the focusing case is above the band-gap frequency and that of the collimated case is below the band-gap frequency. The natural question that may arise is what happens if the frequency traces a broad range of values from below to and above the band-gap of the PC. To investigate this situation, we need first to evaluate the effective index of the PC in the propagation direction (ΓK). The effective index of the PC is calculated in two ways: *n _{eff}* =

*kc*/

*ω*(by PWM) and

where *f* is the air filling factor, *n*
_{1} and *n*
_{2} are the refractive indices of the holes and background, respectively, and *a* is the lattice constant [26, 27]. The values obtained from these calculations are used when a comparison is made between FDTD and Gaussian optics.

We have examined the GRIN PC numerically by FDTD for a range of input beam widths at two representative frequencies above and below the band-gap frequencies (*a*/*λ* = 0.16, *a*/*λ* = 0.24). Figure 6 shows the steady-state field variation as the input beam width *w*
_{0} changes (9.71*a*, 8.92*a*, and 7.93*a*) at *a*/*λ* = 0.16. We can see that the beam converges towards the focal point where the beam amplitude reaches the maximum value *w _{x}* = 4.70

*a*, 4.96

*a*, and 5.48

*a*. The amplitude variation along the

*x*direction at the center of the GRIN PC is displayed in Fig. 7(a). It can be seen that the maximum occurs at the same point which is around 45.5

*a*. (The source is located at 3

*a*distance from the input of the GRIN PC, so it has to be subtracted from the value 48.5

*a*obtained from the graph) We conclude that for the range of

*w*

_{0}listed in Table I, the gradient coefficient

*α*is the same. A larger input beam focuses to a smaller spot-size, and the magnitude variations are consistent with Eq. (12). Similarly, the same calculations are carried out for the frequency at

*a*/

*λ*= 0.24, and Fig. 6 shows the steady-state field variation (

*w*

_{0}takes values 9.38

*a*, 8.53

*a*and 7.65

*a*and

*w*values are 3.25

_{x}*a*, 3.57

*a*, and 4.05

*a*). From Fig. 7(b), the focal point is extracted to be around 40

*a*. [Again, the source is located at 3

*a*distance from the input of the GRIN PC, so it has to be subtracted from the value 43

*a*obtained from Fig. 7(b)] Since the input beam with different widths reaches maximum at the same point, we conclude that

*α*is constant as long as

*w*

_{0}is not too narrow and wide. The calculated values of

*α*from Figs. 7 are 5.75 × 10

^{-4}

*μ*m

^{-1}(

*π*/2

*α*= 45.5

*a*) and 6.54×10

^{-4}

*μ*m

^{-1}(

*π*/2

*α*=40

*a*), respectively. As shown in Fig. 4, the EF curves of a linear PC are expected to behave as a diverging medium for the input beam at these two selected frequencies along ΓK direction (the circular and hexagonal contour shapes are highlighted in the figure). However, the quadratic gradient medium overcomes the diffraction of PC. As a result, the input beams converge toward focal points.

The investigation of the beam profile across *y* cross-section at the focal and input planes provides further insights. Figure 8 indicates the beam profiles at these two planes. The solid and broken lines are for the beams at the input and focal planes, respectively. Table I lists the FWHM values and the conversion spot-size ratios. As the input beam width decreases, the output beam width increases and the spot-size conversion ratio approaches to 1:1, which is the condition for the self-collimation. If the input beam width *w*
_{0} decreases further, then the beam diverges with conversion ratio smaller than one.

To compare the results obtained numerically by FDTD with Gaussian optics, we performed the following calculation. The gradient coefficient *α* is obtained with the values of *w _{x}*,

*w*

_{0}, and

*n*taken from Table I and Fig. 7 for two different frequencies using Eq. (11) (8.64×10

_{x}^{-4}

*μ*m

^{-1}and 8.88×10

^{-4}

*μ*m

^{-1}). Then, with fixed

*α*values, we calculated the output beam width with known input beam width. The calculated values are compared with the ones obtained by FDTD in Table II. The close agreement between the two methods can be seen from the Table II.

Numerically obtained gradient coefficients are smaller than these ones obtained analytically. It is important to note that analytical formulation is carried out for the graded-index slab structure. However, the numerical simulation is performed with a GRIN PC. Even though the air hole refractive indices are modulated, there is another modulated refractive index due to the perforated air holes. The diffractive nature of the PC itself may reduce the power of the focusing. One other reason for the discrepancy is due to the effective index approximation of the GRIN PC. Definitely, such an approximation for a wide range of frequencies contributes to the mismatch. The differences between FDTD simulation (no effective index approximation) and Gaussian optics formulation (with the effective index approximation) at high frequencies are attributed mainly due to the effective index approximation. At high frequencies, Eq. (14) should be replaced by more accurate results obtained from PWM. Despite these factors, the overall characteristics of the GRIN PC can be understood with Gaussian optics and numerical results can be performed to validate the conclusions.

Another question of interest is the cases where the frequency is inside the band-gap (e.g., *a*/*λ* = 0.20) or well above the band-gap (e.g., *a*/*λ* = 0.30). The desirable feature of focusing behavior is lost for the two cases as shown in Fig. 9. For the first, the frequency is inside the band-gap; hence the focusing and propagation are hindered due to PBG effect [Fig. 9(a)]. In Fig. 9(b), on the other hand, the wavelength is almost equal to the lattice constant *a*; hence, the internal diffraction play a role in this case (*a* is comparable with *λ*/*n*
_{eff}). Also the dispersion plot in Fig. 4 supports the observation such that the contour shapes become more complex at higher frequencies. One should notice that operating frequencies at higher bands are not desirable due to the inherit problem of out-of-plane scattering losses (being above the light line).

Finally, to show that the beam can be guided by either self-collimation or oscillatory behavior for a large propagation distance, Fig. 10 is plotted with FWHM values of the input beam taken as 5.5*a* and 7*a*, respectively. It illustrates the field propagation with GRIN PC. It can be seen that the beam width stays constant for the self-collimated case [Fig. 10(a)], and it oscillates for the second case [Fig. 10(b)]. The index gradient produces local modifications of EF curves as shown previously in Fig. 4. Since the beam propagation is along ΓK direction and the dispersion curve is not flat along that direction so we can attribute the observation of the self-collimation at *a*/*λ* = 0.24 case again to graded index medium. Guiding the field for large distances without diffraction is obviously great advantages and oscillatory nature of the beam under appropriate conditions may find interest in the mode matching and coupling/decoupling applications.

## 5. Discussion and conclusion

Below the band gap, graded index medium is responsible for the focusing and self-collimation behaviors as expected. Above the band gap, dispersion effect is expected to play role for the self-collimation. However, as can be seen Fig. 4, the contours are not flat in the propagation direction. So index gradient is the leading mechanism for the self-collimation again. In the band gap, bandgap effect is the leading mechanism. Light tends to bend toward high index medium as it propagates in oscillatory fashion. The input beam width is an initial condition such that it determines the light propagation behavior if it is converging or diverging one at the beginning. However, there is a transition beam width in between the two situations (diverging, converging) and the beam stays collimated if this condition is satisfied. The main reason to study different frequency cases is to underpin the light behavior through all wavelength regimes. The frequency dependency of the focusing and self-collimations comes from the underlying PC structure.

We briefly discuss the feasibility of the fabrication of such a structure operating at 1.55 *μ*m. In the above analysis, either the refractive indices of the air holes or radius are assumed to be modulated. The difficulty of altering refractive indices of holes is obvious at shorter wavelength regime. However, a similar affect can be obtained if the air radius are modulated such that center section has the smallest hole radius, e.g., 0.10*a* in ΓK direction and the radius is increased gradually up to 0.20*a* along ΓM direction with the increment steps of 0.02*a*. Such a fine spatial control may be difficult to fabricate even it is feasible. Instead of PC with *r* = 0.20*a* one can have larger air holes such as 0.30*a* and then the incremental steps can be larger for example 0.03*a*.

In conclusion, we explored for the first time GRIN PCs and investigated the characteristics of it under different conditions. The initial conditions decide if the designed GRIN PC acts as a converging, diverging, or even self-collimated device. Gaussian optics accounts very well for frequencies up to *a*/*λ* = 0.30 at which wavelength is almost equal to the lattice constant. There is deviation from the above categorization for the frequencies which is inside the PBG. The gradient coefficient is frequency dependent, and stays constant for a range of input beam widths. The consistency between numerical result and analytical formulation was presented. The observations obtained for the case study can be generalized for other parameters of triangular-lattice PC. It is known that discrete optical elements to focus or self-collimate the input beam should be integrated for compact optical devices. High level of integration is possible with this approach because there is no need for the discrete alignment. Further optical devices such as input and output couplers can be designed.

## References and links

**1. **J. D. Joannoupoulos, P. R. Villeneuve, and S. Fan, “Photonic crystals: putting a new twist on light,” Nature **386**,143–149 (1997). [CrossRef]

**2. **H. Benisty, J. M. Lourtioz, A. Chelnokov, S. Combrie, and X. Checoury, “Recent advances toward optical devices in semiconductor-based photonic crystals,” Proceedings of the IEEE **94**,997–1023 (2006). [CrossRef]

**3. **B. S. Song, S. Noda, and T. Asano, “Photonic devices based on in-plane hetero photonic crystals,” Science **300**,1537 (2003). [CrossRef] [PubMed]

**4. **O. Painter, R. K. Lee, A. Scherer, A. Yariv, J. D. O’Brien, P. D. Dapkus, and I. Kim, “Two-dimensional Photonic band-gap defect mode Laser,” Science **284**,1819–1821 (1999). [CrossRef] [PubMed]

**5. **B. Saleh and M. Teich, *Fundamentals of photonics* (New York, Wiley, 1991). [CrossRef]

**6. **M. R. Mackenzie and C. Y. Kwok, “Theoretical analysis of integrated collimating waveguide lens,” J. Lightwave Technol. **21**,1046–1052 (2003). [CrossRef]

**7. **M. Zickar, W. Noell, C. Marxer, and N. de Rooij, “MEMS compatible micro-GRIN lenses for fiber to chip coupling of light,” Opt. Express **14**,4237–4249 (2006). [CrossRef] [PubMed]

**8. **A. Mekis, J. C. Chen, I. Kurand, S. Fan, P. R. Villeneuve, and J. D. Joannopolous, “High transmission through sharp bends in photonic crystal waveguides,” Phys. Rev. Lett. **77**,3787–3790 (1996). [CrossRef] [PubMed]

**9. **M. Loncar, T. Doll, J. Vuckovic, and A. Scherer, “Design and fabrication of silicon photonic crystal optical waveguides,” J. Lightwave Technol. **18**,1402–1411 (2000). [CrossRef]

**10. **A. Yariv, Y. Xu, R. K. Lee, and A. Scherer, “Coupled-resonator optical waveguide: a proposal and analysis,” Opt. Lett. **24**,711–713 (1999). [CrossRef]

**11. **M. Bayindir, B. Temelkuran, and E. Ozbay, “Tight-binding description of the coupled defect modes in three-dimensional photonic crystals,” Phys. Rev. Lett. **84**,2140–2143 (2000). [CrossRef] [PubMed]

**12. **H. Kosaka, T. Kawashima, A. Tomita, M. Notomi, T. Tamamura, T. Sato, and S. Kawakami, “Self-collimating phenomena in photonic crystals,” Appl. Phys. Lett. **74**,1212–1214 (1999). [CrossRef]

**13. **X. Yu and S. Fan, “Bends and splitters for self-collimated beams in photonic crystals,” Appl. Phys. Lett. **83**,3251–3253 (2003). [CrossRef]

**14. **D. W. Prather, S. Shi, D. M. Pustai, C. Chen, S. Venkataraman, A. Sharkawy, G. J. Schneider, and J. Murakowski, “Dispersion-based optical routing in photonic crystals,” Opt. Lett. **29**,50–52 (2004). [CrossRef] [PubMed]

**15. **E. Centeno and D. Cassagne, “Graded photonic crystals,” Opt. Lett. **74**,2278–2280 (2005). [CrossRef]

**16. **E. Centeno, D. Cassagne, and J. P. Albert, “Mirage and superbending effect in two-dimensional graded photonic crystals,” Phys. Rev. B **73**,235119–235119 (2006). [CrossRef]

**17. **P. Halevi, A. A. Krokhin, and J. Arriaga, “Photonic crystals as optical components,” Appl. Phys. Lett. **75**,2725–2727 (1999). [CrossRef]

**18. **E. Foca, H. Föll, J. Carstensen, V. V. Sergentu, I. M. Tiginyanu, F. Daschner, and R. Knöchel, “Strongly frequency dependent focusing efficiency of a concave lens based on two-dimensional photonic crystals,” Appl. Phys. Lett. **88**,011102 (1–3) (2006). [CrossRef]

**19. **D. Mori and T. Baba, “Dispersion-controlled optical group delay device by chirped photonic crystals waveguides,” Appl. Phys. Lett. **85**,1101–1103 (2004). [CrossRef]

**20. **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 (1–4) (2006). [CrossRef]

**21. **A. Taflove, *Computational Electrodynamics - The Finite-Difference Time-Domain Method* (Norwood, Massachusetts: Artech House, 2000).

**22. **S. Guo and S. Albin, “Simple plane wave implementation for photonic crystal calculations,” Opt. Express, vol.**11**, pp.167–175, 2003. [CrossRef]

**23. **T. Tamir, *Guided-Wave Optoelectronics* (Springer-Verlag, Berlin, 1990). [CrossRef]

**24. **J. T. Verdeyen, *Laser Electronics* (Prentice Hall, New Jersey, 1995).

**25. **H. Kogelnik, “On the propagation of Gaussian beams of light through lenslike media including those with a loss or gain variation,” Appl. Opt. **4**,1562–1569 (1965). [CrossRef]

**26. **S. M. Rytov, “Electromagnetic properties of a finely stratified medium,” Sov. Phys. JETP, **2**,466–475 (1956).

**27. **P. Lalanne, “Effective medium theory applied to photonic crystals composed of cubic or square cylinders,” Appl. Opt. **27**,5369–5380 (1996). [CrossRef]