Analyzing the propagating waves in the two-dimensional photonic crystal by the decoupled internal-field expansion method

decoupled internal-field expansion method

Ting-Hang Pei and Yang-Tung Huang

Citation: AIP Advances 2, 012188 (2012); doi: 10.1063/1.3700236

View online: http://dx.doi.org/10.1063/1.3700236

View Table of Contents: http://scitation.aip.org/content/aip/journal/adva/2/1?ver=pdfcov Published by the AIP Publishing

Articles you may be interested in

Influence of unit cell rotated on the focusing in a two-dimensional photonic-crystal-based flat lens J. Appl. Phys. 101, 123112 (2007); 10.1063/1.2745374

Polarization beam splitter based on a two-dimensional photonic crystal of pillar type Appl. Phys. Lett. 89, 171115 (2006); 10.1063/1.2360201

Line-defect-induced bending and splitting of self-collimated beams in two-dimensional photonic crystals Appl. Phys. Lett. 87, 181106 (2005); 10.1063/1.2112186

Time-domain measurement of picosecond light-pulse propagation in a two-dimensional photonic crystal-slab waveguide

Appl. Phys. Lett. 84, 4690 (2004); 10.1063/1.1760224

Rapid photon flux switching in two-dimensional photonic crystals Appl. Phys. Lett. 84, 1817 (2004); 10.1063/1.1667593

Analyzing the propagating waves in the two-dimensional

photonic crystal by the decoupled internal-field expansion

method

Ting-Hang Pei1,a_{and Yang-Tung Huang}2

1_{Microelectronics and Information Systems Research Center, National Chiao Tung}

University, 1001 Ta-Hsueh Road, 30010, Hsinchu, Taiwan

2_{Department of Electronics Engineering and Institute of Electronics, National Chiao Tung}

University, 1001 Ta-Hsueh Road, 30010, Hsinchu, Taiwan

(Received 6 December 2011; accepted 9 March 2012; published online 27 March 2012)

We propose the decoupled internal-field expansion (IFE) method to discuss refrac-tions in the photonic crystal (PhC). This method decouples the full wave in the PhC and classifies them into two categories including the forward-propagating and backward-propagating waves denoted by the index n. A triangular-PhC case is demon-strated and both positive and negative refractions are discussed by this method. The incident angle of 10◦results in the positive refracted wave with the refracted angle about 8◦, which approximately corresponds to the forward wave of n=0 order. The negative refracted waves, which exist in the left and right edge regions, propagate almost parallel to the interfaces between the PhC and outside media. Meanwhile, due to the interaction between the negative refracted wave and the nearest few rows of air cylinders, the reflected wave and another weaker negative refracted wave are created. Finally, the weaker negative refracted waves from both edge regions inter-fere with each other in the middle region. It is found out that the negative refracted waves in edge regions as well as the interfered wave in the middle region can be con-structed by two n=-1 forward and backward waves. On the other hand, the positive refracted wave is composed of the n=0 forward wave dressed other n=0 forward waves, the propagation angle is affected by these dressed waves, especially near the edge region. Finally, another case proves this point of view explicitly. Copyright 2012 Author(s). This article is distributed under a Creative Commons Attribution 3.0 Unported License. [http://dx.doi.org/10.1063/1.3700236]

I. INTRODUCTION

It is well known that the photonic crystal (PhC) formed with dielectric periodic structures, and the Bloch theorem can describe the electromagnetic characteristics very well, which is analogous to the case of the electron in solids. The photonic band structure (PBS) of the PhC reveals lots of information on the propagation of light. One of them is the light confinement due to the existence of the photonic band gaps (PBGs).1,2_{It induces some applications such as photonic crystal waveguides,}

which can bend light extreme efficiently around sharp corners and have been successfully developed in two-dimensional and three-dimensional structures.3–6

Recently, two popularly refractive phenomena, the so-called the superprism effect and the negative refraction in the PhC, have been investigated and discussed.7–14All of them treat PhCs as homogeneous media with negative refractive indices, in which the propagation directions of EM waves can be determined from the equifrequency surface (EFS). The normal direction of the EFS is the direction of the group velocity as well as the propagation direction of light in the PhC, so the EFS can give the relation between the incident and the refracted angles without considering the interface effect.15_{No matter the compositions of the PhC are, Snell’s law can be used to describe single-beam}

a_{thp3000.ee88g@nctu.edu.tw}

refraction between the PhC and other homogeneous medium.17_{Although the EFS can predict a lot}

of variant optical performances such as the superprism effect and the negative refraction, however, the PhC is always finite in practice. After all, the EFS represents the characteristic of an infinite PhC, not a finite one. Furthermore, the existence of interfaces may affect optical performances which are different from the predictions of the infinite PhC. It causes us to design some optical devices carefully when the PhC is used.

The affection of the interface on the propagation direction of light in the PhC has been proved. It is the fact that the symmetry is broken and the multiple scattering no longer satisfies the periodic condition at the interface. So the full wave inside a finite PhC cannot be only composed of Bloch waves. Felbacq et al.16_{shows that the wave propagating through the PhC is dressed by evanescent}

waves near the interface of the PhC with finite size. They use the concept that photons strongly scattered by the boundary of the PhC resulting in the creation of evanescent waves. Finally, the refracted angle has to be fixed. They classify the interface effect as a generalized version of the Goos-Hanchen effect.

On the other hand, the phase matching condition, which is widely used to study 1D diffraction gratings, has also been used to discuss refraction at an interface between a homogeneous medium and a 2D PhC very well.13,14,17,18_{Under this condition, the tangent component of the incident wave}

vector has to equal the sum of the tangent components of refracted waves. It can be extendedly applied to the wave vectors located in other primitive cells, not only within first Brillouin zone (BZ). Furthermore, anomalous refractive phenomena in the PhC can also be explained under this condition with the wave vector diagram formalism.13,14,17_{Usually, the refracted wave lying in the first BZ is a}

zeroth-order wave, and others belong to higher Brillouin zones are the higher-order waves. Either the zeroth-order wave or higher-order ones have been investigated by the finite-difference time-domain (FDTD) simulations.13,14,17,18_{However, the waves in the PhC include not only propagation waves,}

but also evanescent waves. The later is hard to be filtered out from the FDTD method. As mentioned before, the effective refractive index of the finite PhC is different from that of the infinite PhC due to the excitation of evanescent waves at the interface. If we want to understand how evanescent waves affect the propagation waves and the field distributions of them, we have to use another method other than the FDTD one.

As we know, the physics of the refraction is due to the coupling between the incident EM wave and one or more propagation waves in the PhC.13,14,17,18_{Sakoda has successively calculated}

the transmission and the Bragg reflection spectra using the internal-field expansion (IFE) method and mentioned existences of the uncoupled modes.19–21_{He also pointed out that the strength of the}

coupling between the outside waves and the internal field at the interface substantially modifies the transmission.19 _{This method has been widely used in studying the PhCs with graphite structure}22

and the vanishing transmission due to symmetry.23–33_{But this method shows the field distribution}

as sum of all expanded series which is also equal to sum of all propagating waves displayed by the FDTD method. It cannot exhibit right information inside the PhC if all series are combined together. However, after rearranging some equations, we find out that this method can discuss each propagating wave shown by the FDTD method and how evanescent waves affect propagating ones in the PhC as Felbacq et al. revealed.16_{We call this new way as the decoupled internal-field expansion}

method.

In this paper, both the FDTD method and the decoupled IFE method are used to discuss refractions in the PhC. By the FDTD method, we can investigate the evolution process of propagating waves inside the PhC. When we discuss the refraction by the decoupled IFE method, each propagating wave can be reconstructed by summing of some expanded series. We also discuss the effect of interfaces on propagating waves due to evanescent ones. In addition, the effect of the width in the edge region is also brief discussed.

II. THEORY AND NUMERICAL METHOD A. Fields inside photonic crystal

We consider a two-dimensional PhC composed of triangular array of air cylinders with radius of r in a dielectric background. The dielectric constants of the cylinders and the background are

x

a

w repetitions

k

k

k

k

2r

Air

Air

Supercell

Transmitted

Region

Incident

Region

y

d

d

k

k

k

εaandεb, respectively. The infinitely extended direction of air holes is parallel to the z-axis. The PhC is infinitely extended in the x-direction and has finite width in the y-direction. Therefore, two interfaces between the dielectric and the PhC are located at x=0 and x=L, respectively, where L is the total width of the PhC in the y-direction. The region from the interface to the nearest cylinder is called the edge region. There are two edge regions in this structure. They are the left and right edge regions with widths of d1 and d2, respectively. The region between two edge regions is called the middle region. In order to describe this structure, a1and a2are defined as two lattice periods parallel and perpendicular to the interface, respectively. Then the total width is L= (2r+d1+d2)+(w-0.5)a2, where w is the number of repetitions in the y-direction. The configuration of the PhC is shown in Fig.1, and the first BZ is also shown in the up-right corner.

Next, let us consider an E-polarized (electric field parallel to the z-axis) plane wave incident upon the left interface along theM direction. The 2D wave vector of the incident wave is denoted byki = (kisinθ, kicosθ) = (ki,x,, ki,y), whereθ is the incident angle and ki = √εiω/c. In the representation of ki,εiis the dielectric constant of the incident region,ω is the angular frequency of the incident wave, and c is the light velocity in vacuum. The incident wave would excite some propagating modes in the PhC and cause different refractions. We follow the IFE method19–21_{and extend its function to}

discuss refractions in the PhC. This method is based on internal fields expanded in Fourier series, and the electromagnetic (EM) wave in the incident region is represented as a Rayleigh expansion constituting the incident plane wave and the reflected Bragg waves. In the transmitted region, the EM wave is composed of transmitted Bragg waves. The reflected and transmitted Bragg waves of order n are represented as space harmonics with the wave vector

knr_,x = k n t_,x = k n x = ki_,x + nG, (1) where kn

r_,x and ktn_,x are the components of wave vector parallel to the interface for the reflected and transmitted Bragg waves, respectively, and n is an integer. G=2π/a1 is the reciprocal lattice

vector corresponding to the surface periodicity a1. Each component of Eq.(1)is called the nth order phase matching condition. The wave vector components normal to the interface for the reflected and transmitted Bragg waves of order n are

krn_,y= ⎧ ⎨ ⎩ −k2 i −  kn x 2 if ki ≥kxn −ikn x 2 − k2 i otherwise, (2) ktn,y= ⎧ ⎨ ⎩ −kt2−  kn x 2 if kt ≥kxn −ikn x 2_{− k2} t otherwise, (3)

Here, kt = √εtω/c and εtis the dielectric constant in the transmitted region.

The electric fields in the incident region and the transmitted region are respectively given by Ei z(x, y) = E0exp(i  ki·r )+ ∞ n_=−∞ Rnexp(i  knr ·  r ), (4) Et z(x, y) = ∞ n_=−∞ Tnexp(i  knt · _{r ).} (5) The electric field in the PhC is expanded as

EPC(x, y)= ∞ n=−∞ _y LTn+  1−y L  (δn0E0+ Rn) exp(i k_xnx)+ ∞ n=−∞ ∞ m₌₁ Anmexp(i knxx) sin mπ L y, (6) which satisfies the following equation:

 1 ε(x, y)  _∂2 ∂x2 + ∂2 ∂y2  +ω2 c2  EPC(x, y) = 0. (7)

The inverse ofε(x, y) is also expanded in Fourier series within a symmetrical supercell region:19–21

1 ε(x, y) = ∞ n_=−∞ ∞ m=1 κnmexp(i Gnx+ kmy), (8)

where km=mπ/L. Following the same calculation processes19–21and consider boundary conditions for E-polarized wave, we can determine the unknown coefficients, Anm, Rn, and Tnwhere Rnand Tn are relate to the reflection and transmission spectra, respectively. In practical calculation, the Fourier series of the ε(x, y) are truncated and expanded into (2N+1)M terms. Then the total number of the unknown coefficients is (2N+1)(M+2). From Eqs.(4)and(5)and boundary conditions, we also obtain (2N+1)(M+2) independent equations. Solving these independent equations, all coefficients can be obtained.

Due to two parallel interfaces, multiple scatterings take place in the finite PhC. It leads that electric fields can be decomposed into two sets. One is related to exp(ikmy), and the other is related to exp(-ikmy). The former is classified as the forward-propagating wave (+y-direction) and the latter the backward-propagating wave (-y-direction).16 For convenience, they are shortly called the forward wave and the backward wave in the following discussion. About these two kinds of waves, Jiang et al.18 pointed out that the number of forward waves is equal to the number of backward waves. Simply expanding sinkmy by (exp(ikmy)-exp(ikmy))/2i in Eq.(6), the Ez-fields of the forward and backward waves are then given by

Ezn+ = 1 i ∞ m=−∞  Anm 2 +  (−1)m−1_T n+ Rn+ δn0E0  mπ  × expiknxx+ kmy  , (9)

Ezn− = − 1 i ∞ m_=−∞  Anm 2 +  (−1)m_−1T n+ Rn+ δn0E0  mπ  × expikxnx− kmy  . (10)

Furthermore, the corresponding magnetic fields are Hxn+ = −i ωμ0 ∞ m=−∞  Anm 2 +  (−1)m−1_T n+ Rn+ δn0E0  mπ  ×mπ L  expikxnx+ kmy  , (11) Hyn+ = i ωμ0 ∞ m_=−∞  Anm 2 +  (−1)m_−1T n+ Rn+ δn0E0  mπ  × kn xexp  ikxnx+ kmy  , (12) Hxn− = − 1 i ∞ m_=−∞  Anm 2 +  (−1)m₋₁ Tn+ Rn+ δn0E0  mπ  ×mπ L  expiknxx− kmy  , (13) H_yn− = −i ωμ0 ∞ m=−∞  Anm 2 +  (−1)m−1_T n+ Rn+ δn0E0  mπ  × kn xexp  ik_xnx− kmy  , (14)

where μ0 is the magnetic permeability in vacuum. Once these coefficients Anm, Rn, and Tn are known, the forward and backward waves can be obtained by substituting these coefficients into Eqs.(9)–(14). We call this extended method as the decoupled Internal-field expansion method.

B. Group velocity and Poynting vector

The propagation direction inside the PhC is commonly determined by using the group velocity vectorvg = ∇

kω.

34,35_{It has been proved that}3,34



ve =



vg= ∇

kω, (15)

whereve is the energy velocity in the periodic medium and defined as  ve = 1 S   P kd  r_// 1 S  U kd  r_// ≡ P k U k , (16)

where S is the area of a unit cell,P

kis the time-averaged Poynting vector, and Uk is time-averaged energy density. The unit cell, which is the rectangular region closed by dashed line, is shown in Fig.1. Eq. (15) means that we can obtain the propagation direction inside the PhC by calculating the averaged energy velocity in a unit cell including the area occupied by the scatters.

The time-averaged Poynting vector at spatial point is given by

 P k = 1 2Re  E k(  r )×H∗ k(  r ). (17)

We drop the words “time-averaged” and call P

k the Poynting vector for convenience. Inside the PhC, the x and y components ofP

k can be written as P k,x = −1 2 Re ∞ n=−∞  E+znHyn+∗+ E−znHyn−∗  , (18)

P k_,y= 1 2Re ∞ n_=−∞  E+znHxn+∗+ E−znHxn−∗  . (19)

Eqs.(18)and(19)can be further used to discuss the contribution of the nth-order wave or the sum of several waves. Finally, the propagation angleθris defined as follows:

tanθr = P k,x P k,y . (20)

III. ANALYSIS OF THE REFRACTIVE PHENOMENA A. The FDTD Simulations

The Bragg transmission and reflection spectra of the structure in Fig.1have been calculated by Sakoda.19–21The refraction in this plane-parallel slab is not like that in the semi-infinite PhC with a single interface. No matter how far the distance between two interfaces is, multiple reflections always take place which lead to reflected waves interfere with the original propagating wave and other reflected waves. However, due to the finite width of the incident beam like the Gaussian beam in the FDTD simulation, the interference may not overlap in all space even internal reflections successively produce. It tells us that we have to design practical devices carefully because propagations inside the finite PhC may be deviated the realized situation of the plane wave incidence.36 _{So in this section,}

we first use the FDTD method to investigate refractions in the finite PhC and explain some special phenomena physically, and then we discuss the compositions and the propagation direction of each propagating wave by the decoupled IFE method. Finally, the results of the decoupled IFE method are compared with those of FDTD simulations.

We follow the designed parameters in Ref.19to investigate refractions and field distributions inside the PhC. The radius of each air hole is r=62.5 μm and the lattice constant a2= a1/

√ 3 is 170μm, which implies r/a2=0.367647. The dielectric constants are εa=1.0, εb=2.1, εi=1.0, and εt=1.0. The frequency f of the incident beam is 0.375 (c/a2) in the first photonic band where the gradient of the EFS points outward. In our discussions, both interfaces are extended along theM direction. According to EFSs of the first photonic band, the relation between the incident angle and the refracted angle at 0.375 (c/a2) is shown in Fig.2. Total 5041 plane waves are used in the EFS calculations and the error is believed below 1%. In this case, two relation carves exit. One describes the positive refraction, and the other represents the negative one. The positive refraction takes place from the incident angle 0◦and the negative one can be observed above 16◦. Although the negative refraction starts from 16◦, the refracted wave which is attenuated along the interface can be still found out below the incident angle 16◦. The incident angle is chosen as 10◦for the following discussions.

Then the FDTD method is used to investigate refractions in the finite PhC. The number of repetitions w along the y-direction is 29 and both widths of the left and right edge regions are equal to 2a2. A E-polarized Gaussian beam is incident from air into the PhC with the incident angle 10◦. Two spatial field distributions at different times are shown in Fig.3. One at earlier time is drawn in Fig.3(a), and the other at later time is shown in Fig.3(b). The incident regions are at the left in both figures. In Fig.3(a), the incident wave excites two main waves in the PhC. One much larger wave propagates through the PhC in the manner of positive refraction with the refracted angle about 8◦. Another small one exists in the left edge region and then propagates along the interface. This wave seems to be confined in the left region and behaves as a negative refracted wave which is not predicted in Fig.2. Further investigations from Figs.3(a)and3(b)show that some waves leak out the interface and radiate back into the incident region like the reflection wave. By investigating the field distribution at the early time in Fig.3(a), it shows that the negative refracted wave brings about these leaked waves. Although the left region is like a waveguide sandwiched between air and the PhC, however, energy is not totally confined. The magnitude of the Poynting vector shown in Fig.4, which is calculated along the interface in the left region and averaged in the area like the

FIG. 2. The relation between the incident angle and the refracted angle at 0.375 (c/a2).

supercell in Fig. 1, displays that the energy flow decreases when the wave propagates far away from the incident place. So the negative refracted wave in the left region is unstable and decays as the distance increases. After propagating through about 15a1 and 180a1, the Poynting vector decreases to 90% and e-1_{of the original strength, respectively. The cause is the interaction between} the negative refracted wave and the nearest few air-hole layers. The evidence can be seen in Fig.3

where the reflected wave produces as the negative refracted wave propagates forwardly. Besides, this interaction also successively induces the second kind of negative refracted waves which enters the middle region and propagates toward the right edge region as the left arrow denotes in Fig.3(a). This negative refracted wave is much weaker than the one in the left edge region. Its propagation direction forms an included angle about 10◦with the x-direction. Eventually, this wave touches the right edge region and a certain part reflects to the middle region. By the way, it tells us that this structure could be severed as a well function of the waveguide if the length in the x-direction is only few dozen of a1.

Next, the wave in the right edge region is discussed. Because the compositions of the PhC are all lossless in our discussions, the eigenfrequencyω is real which implies the dispersion relation possessing the time-reversal symmetry.37 According to this symmetry, the process of the incident wave exciting the positive and negative refracted waves can be reversed. When the positive refracted wave leaves the middle region into the right edge region, another wave and its propagation direction like the negative one in the left region can be generated. In Fig.3(a), the wave in the right region takes place later than the one in the left region due to causality. Like the phenomenon in the left region, the wave in the right region also interacts with the nearest few air-hole layers and then leaks out some waves into the transmitted and middle regions. So this wave is also unstable as the one in the left region. The induced wave enters the middle region and propagates in the direction as the arrow denotes. Then the induced waves from the left and right regions interfere each other in the middle region. The circled part in Fig.3(a)shows that the interfered wave is in the process of forming. Finally, the interfered wave establishes a special field distribution as the circled part shows

d1 d2 Middle Region (a) 0 30a1 60a1 90a1 120a1 150a1 0 30a1 60a1 90a1 120a1 150a1 d1 Middle Region d2 (b)

FIG. 3. The field distributions of (a) the early time and (b) the later time in the case of d1=d2=2a2,εa=1.0, εb=2.1, εi=1.0,

FIG. 4. The magnitude of the Poynting vector along the interface in the left edge region.

in Fig.3(b). It can be seen that a lot of similar field distributions repeat along the x-direction in the middle region. In this case, special phenomenon such as the negative refraction and the interfered wave can be investigated even the incident angle is less than 16◦. So the EFSs cannot always predict some propagation phenomena in the finite PhC due to the broken symmetry. These phenomena were very rarely mentioned in past references.

At the end of this section, let us shortly discuss the effect of the edge region. It is found out that the whole field profile can be varied if the widths of both edge regions are changed. In Fig.5, both d1and d2are changed to 7a2. As mentioned before, the function of each edge region is like a waveguide confined in between two mediums if the PhC is treated as an effective medium. As the width of the edge region is changed, the wave in it is also altered. It is just like the different mode in the waveguide if the width of the waveguide is broadened or narrowed. Besides, the propagation direction is also different from the case of d1=d2=2a2. By comparing Fig.5with both Figs.3(a) and3(b), the propagation direction in the edge region in Fig.5is deviated the x-direction more than that in Fig.3. The interfered wave in the middle region, which is not regularly repeated along the x-direction, is also different from that in Fig.3. It is the reason that the width of the edge region changes the propagation direction of the wave in the edge region as well as that of the induced wave in the middle region. In conclusion, some part of the field distribution and the radiation outside the PhC can also be affected by the width of the edge region. If both widths of the left and right regions are asymmetry, the field distribution would be much more complicated than the symmetry case.

B. The Decoupled IFE Calculations

Although the FDTD simulations can give real-time information, some effects like the evanescent waves near the edge are hard to know through the FDTD method. Furthermore, the grating-like structure such as the PhC includes some special propagating waves. The optical performance in each of them should be analyzed individually in order to understand more characteristics of the PhC. By deeply comprehending the IFE method, it is possible to solve the above problems through this method. This method contains the total wave in the PhC including the forward and backward waves which cannot easily be pointed out by the FDTD simulations. It has been introduced in detail in Sec.II.

Middle Region d1 d2 0 30a1 60a1 90a1 120a1 150a1

By using Eqs.(9)–(14), the total wave inside the PhC can be decoupled into infinite sub-fields. One advantage of this method is to reveal the relation between each sub-field and the width of the structure in the y-direction. It can tell us how each sub-field varies with the width L. In order to verify our extended method, we repeat the same case in section A. The frequency is still 0.375 (c/a2), the incident angle is 10◦, and other parameters remain. The expanded series in Eqs.(9)–(14) are truncated at -5n5 and 1m260 leading to total 2882 unknown coefficients. The choices of n and m are enough for convergence.17–19_{After solving the 2882×2882 coefficient matrix, each}

sub-field can be obtained by substituting Anm, Rn, and Tninto Eqs.(9)–(14).

First, three propagating waves of zeroth order (n=0), forwardly minus first order (n=-1), and backwardly minus first order (n=-1) are demonstrated in Figs.6(a)–6(c). Both sides in Fig.6are two interfaces between the PhC and air, so the field distributions in the incident and transmitted regions are not shown. In Fig.6(a), the n=0 forward wave like a plane wave propagates in an effective

medium with an tilted angle of 8◦ related to the y-direction, which is calculated by Eq.(20). This wave is approximated to the positive refracted waves about between 20a1to 45a1in the x-direction in Figs.3(a)and3(b). The color contrast in Fig.3is not so sensitive like that in Fig.6(a)because we want to stress the importance of the interfered waves in the middle region in Fig.3. The n=0 backward wave has the same field profile by rotating the n=0 forward wave counterclockwise 180◦, but the Poynting vector of the backward one is only about 5% as much as that of the forward one. The field distributions of the forward and backward n=-1 waves are shown in Figs.6(b)and6(c), respectively, where the most part of them exist in the left and right edge regions and very few in the middle region. They are almost the same if we rotate one of them 180◦counterclockwise.

Next, some former figures in Fig.6are combined to form new figures. Then these new figures can be used to represent several parts of field distributions in the FDTD simulations. The first new figure is shown in Fig.6(d), which is produced by combining Figs.6(a)–6(c)and the n=0 backward

wave. This figure corresponds to the field distribution about between 40a1to 45a1in the x-direction of Fig.3(b). If we enlarge the corresponding region in Fig.3(b), it is very similar to Fig.6(d)whether in the edge region or the middle region, especially the field patterns in two edge regions. It means that if the backward waves are filtered out from Fig.7(d), the residual wave can represent the positive refracted wave about between 20a1to 45a1precisely.

The second new figure is shown in Fig.6(e), which is created by summing Figs.6(b)and6(c). In this figure, the most part of fields exist in both edge regions, which are the same as those in Fig.3(b). By using Eqs.(18)–(20), the propagation directions in both edge regions are almost parallel to the x-direction. Besides, it can be investigated that the field distribution repeating in the middle region is also alike that in Fig.3(b)as mentioned before. Summary, the results of propagation directions and the field distribution shows that the combination of the n=-1 forward and backward waves gives both waves in two edge regions of the FDTD result in Fig.3(b).

C. Propagation Angles

Since the previous section shows that the decoupled IFE method can explain the FDTD simu-lations and construct the field distributions of different propagating waves, the results permit us to calculate the propagation angles by Eqs.(18)–(20). The Poynting vector as well as the propagation angle are calculated and averaged in a supercell as shown in Fig.1. Starting at the first air-hole layer close to the left edge region, the supercell is numbered along the y-direction. Finally, there are total 29 supercells in which the last one includes parts of the right edge region. We only consider and discuss the results of the first 28 supercells due to the above reason.

In Fig.7(a), four results about the forward and backward waves of n=0 and n=-1 orders are shown. The propagation angles of the n=0 forward and backward waves in each supercell are kept at almost the same values even close to the edge regions. On the contrary, the propagation angles have large variances near edge regions for both n=-1 forward and backward waves. In the case of the n=-1 forward wave, the propagation angle is close to -90◦_{at the 1st and 28th supercells, and the} absolute values decrease from the edge supercell toward the central one. The explicitly decreasing regions cover the left and right outermost 9 supercells. At the 16th supercell, the propagation angle tends to be a stable value of about -29◦. But the situation for the backward wave has something

(a)

(b)

(c)

(d)

(e)

FIG. 6. The field distributions of (a) the n=0 forward wave, (b) the n=-1 forward wave, (c) the n=-1 backward wave, (d) the combined wave of n=0 and n=-1 of the forward and backward waves, and (e) the combined wave of the n=-1 forward and backward waves.

(a)

(b)

0 10 20 30 The Number of Supercell

10 0 -20 -40 -60 -80 0 10 20 30 The Number of Supercell

0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 (c) 0 10 20 30 The Number of Supercell

0 10

-10

FIG. 7. (a) The propagation angles and (b) the Poynting vectors of n=0 and n=-1 of the forward and backward waves, and (c) the propagation angles of the n=0 forward wave, the total forward wave, and the full wave in each supercell.

FIG. 8. The field distribution in the case of the frequency 0.52 (c/a2) and the incident angle 10◦. Other parameters are the

same as those in Fig.3.

different with the forward one. Its propagation angle is only -75◦less than that of the forward wave. The decreasing regions only cover the left and right outermost 4 supercells. So the backward wave reaches a stable condition faster than the forward one. The stably absolute value of the backward wave is about -44◦large than that of the forward one. From Fig.7(b), the Poynting vectors of the n=-1 forward and backward waves in each supercell are almost the same. For two n=0 waves, the Poynting vector of the forward one is twenty times as much as that of the backward one. So in the interfered region including two n=0 waves, the propagation angle is mainly determined by the forward wave. Besides, the forward n=0 wave carries the transmitted information and the backward one the reflected information. The ratios of Poynting vectors between two n=0 waves also means that the transmission and reflection in this case are about 95% and 5%, respectively, which are very close to the FDTD calculations.

Finally, the propagation angles of the n=0 forward wave, the total forward waves, and the full wave are shown in Fig.7(c), where three kinds of waves exists mainly about between 20a1 to 45a1in the x-direction of Fig.3(b). The full wave includes both forward and back waves, in which the main part displays as Fig.6(d). It can be found out that the propagation angles of these three waves are very close to each other in all supercells. Comparing the n=0 forward wave with the total forward waves, all other n=0 forward waves contribute to the two supercells nearest to the edge regions, which result in the propagation angles a little reduction. If the interfaces disappear, the propagation angle of the n=0 forward wave is the same as that calculated from EFSs of the infinite PhC. So the broken symmetry creates some waves mainly near the interfaces. The result is the same as Felbacq shown,16 _{in which the evanescent waves cause the propagation angles of the}

Bloch waves to decrease in the spaces near the interfaces between the periodic-cylinder region and the homogeneous media. In our discussion, all other n=0 and n=-1 waves construct the evanescent waves in the PhC. In Fig.8, we show another case in which the frequency is 0.52 (c/a2), the incident angle is 10◦, and all other parameters are the same as those in Fig.3. In this case, only one refracted beam is investigated. When the n=0 forward and backward waves are filtered from the full wave, the residual field distribution mainly exist at the places near two edge regions, which represents the evanescent wave that affects the propagation angles of the n=0 forward and backward waves at these places. By using the decoupled IFE method, the evanescent wave can be clearly constructed which is hard to be extracted from the FDTD simulation. Another information in Fig.7(c)is that the total backward waves cause the propagation angles of the full wave a little increase at the one-fourth and three-fourth middle region and their neighboring places.

IV. CONCLUSION

In this paper, the decoupled IFE method can give some information which cannot be found out easily through the FDTD method. In one analyzed case, the positive refracted wave shown in the FDTD calculation can be approximately constructed by summing of the n=0 and n=-1 forward

waves. Its propagation angle is the same as the n=0 forward wave when the supercells are away from the edge regions. The deep place in the finite PhC is close to the situation of the infinite PhC, where the propagation angle of the finite PhC fits the EFS result. The difference between the finite and infinite PhCs in the propagation phenomena is arising from the broken symmetry due to the interfaces. We prove that the propagation angles of the supercells near the edge regions have a little decrease due to all the other n=0 forward waves, which form evanescent waves.

Besides proposing the decoupled IFE method, some interesting phenomena are also discovered in this paper. We find out that the negative refracted waves exist, which is not predicted by the EFSs but can be exhibited by the FDTD method as well as the decoupled IFE method. The negative refraction should start at the incident angle 16◦, not 10◦in our case. However, the negative refracted waves, which mainly exist in the left and right edge regions and propagate almost parallel to the interfaces between the PhC and outside media, are attenuated ones. They are unstable until the incident angle is 16◦. Due to the interaction between the negative refracted wave and the nearest few rows of air cylinders, the reflected wave and another weaker negative refracted wave with a small refracted angle are created. The weaker negative refracted waves produced from the left and right edge regions together form the interfered wave in the middle region. It is found out that the negative refracted waves in both edge regions as well as the interfered wave in the middle region can be approximately constructed by two n=-1 forward and backward waves.

In conclusion, the decoupled IFE method responds to the FDTD simulations and analyzes the compositions of each propagating wave. It compensates the insufficient predictions of EFSs due to finite size or broken symmetry. Furthermore, the field distributions of evanescent waves inside the PhC can also be given and explain the deviation of the propagation angle near interfaces. So this method has the ability to exhibit important information when discussing the refractions in the finite PhC.

1_{S. John,Phys. Rev. Lett.58, 2486 (1987).}

2_{A. Z. Genack and N. Garcia,Phys. Rev. Lett.66, 2064 (1991).} 3_{A. Yariv, Y. Xu, R. K. Lee, and A. Scherer,Opt. Lett.24, 711 (1999).} 4_{M. Bayindir, B. Temelkuran, and E. Ozbay,Phys. Rev. Lett.84, 2140 (2000).}

5_{Steven G. Johnson, Pierre R. Villeneuve, Shanhui Fan, and J. D. Joannopoulos,Phys. Rev. B62, 8212 (2000).}

6_{M. Bayindir, E. Ozbay, B. Temelkuran, M. M. Sigalas, C. M. Soukoulis, R. Biswas, and K. M. Ho,Phys. Rev. B63, 081107}

(2001).

7_{H. Kosaka, T. Kawashima, A. Tomita, M. Notomi, T. Sato, and S. Kawakami,Phys. Rev. B58, R10096 (1998).} 8_{S. Foteinopoulou and C. M. Soukoulis,Phys. Rev. B67, 235107 (2003).}

9_{Natalia Malkova, David A. Scrymgeour, and Venkatraman Gopalan,Phys. Rev. B72, 045144 (2005).} 10_{Suxia Yang, Tao Xu, and Harry Ruda,Phys. Rev. B72, 075128 (2005).}

11_{R. Moussa, S. Foteinopoulou, Lei Zhang, G. Tuttle, K. Guven, E. Ozbay, and C. M. Soukoulis,Phys. Rev. B71, 085106}

(2005).

12_{Irfan Bulu, Humeyra Caglayan, and Ekmel Ozbay,Phys. Rev. B72, 04514 (2005).} 13_{S. Foteinopoulou and C. M. Soukoulis,Phys. Rev. B72, 165112 (2005).} 14_{Alejandro Martinez and Javier Marti,Phys. Rev. B71, 235115 (2005).} 15_{M. Notomi,Phys. Rev. B62, 10696 (2000).}

16_{Didier Felbacq and Rafik Smaali,Phys. Rev. Lett.92, 193902-1 (2004).} 17_{Israel De Leon and Filippus S. Roux,Phys. Rev. B71, 235105 (2005).} 18_{Wei Jiang, Ray T. Chen, and Xuejun Lu,Phys. Rev. B71, 245115 (2005).} 19_{Kazuaki Sakoda,Phys. Rev. B51, 4672 (1995).}

20_{Kazuaki Sakoda,Phys. Rev. B52, 8992 (1995).}

21_{Kazuaki Sakoda, Optical Properties of Photonic Crystals (Springer, Berlin, 2001).} 22_{D. Cassagne, C. Jouanin, and D. Bertho,Appl. Phys. Lett.70, 289 (1997).}

23_{S. Rowson, A. Chelnokov, C. Cuisin, and J.-M. Lourtioz,J. Opt. A: Pure Appl. Opt.1, 483 (1999).} 24_{S. Rowson, A. Chelnokov, and J.-M. Lourtioz,J. Lightwave Tech.17, 1989 (1999).}

25_{H. Benisty, C. Weisbuch, D. Labilloy, M. Rattier, C. J. M. Smith, T. F. Krauss, R. M. de la Rue, R. Houdre, U. Oesterle,}

and C. Jouanin,J. Lightwave Tech.17, 2063 (1999).

26_{J. Schilling, R. B. Wehrspohn, A. Birner, F. M¨uller, R. Hillebrand, U. G¨osele, S. W. Leonard, J. P. Mondia, F. Genereux,}

H. M. van Driel, P. Kramper, V. Sandoghdar, and K. Busch,J. Opt. A: Pure Appl. Opt.3, S121 (2001). 27_{Ying Xu, Hong-Bo Sun, Jia-Yu Ye, Shigeki Matsuo, and Hiroaki Misawa,J. Opt. Soc. Am. B18, 1084 (2001).} 28_{Anne-Laure Fehrembach and Anne Sentenac,J. Opt. Soc. Am. A20, 481 (2003).}

29_{Anne-Laure Fehrembach and Anne Sentenac,Appl. Phys. Lett.86, 121105 (2005).}

30_{Onur Kilic, Michel Digonnet, Gordon Kino, and Olav Solgaard,Opt. Express16, 13090 (2008).} 31_{Ken-Ming Lin and G. Y. Guo,J. Opt. Soc. Am. B25, C75 (2008).}

32_{Ilaria Gallina, Marco Pisco, Armando Ricciardi, Stefania Campopiano, Giuseppe Castaldi, Andrea Cusano, and}

Vincenzo Galdi,Opt. Express17, 19586 (2009).

33_{Massimo Moccia, Marco Pisco, Antonello Cutolo, Vincenzo Galdi, Pierantonio Bevilacqua, and Andrea Cusano,Opt.} Express20, 18842 (2011).

34_{A. Yariv and P. Yeh, Optical Waves in Crystals (Jonh Wiley & Sons, 1984).} 35_{Chao-Hsien Kuo and Zhen Ye,Phys. Rev. E70, 046617 (2004).}

36_{Z. Y. Li and K. M. Ho,Phys. Rev. B68, 155101 (2003).}