## Abstract

Using the three-dimensional (3D) finite-difference time-domain (FDTD) method, we have investigated in detail the optical properties of a two-dimensional (2D) photonic crystal (PC) surface-emitting laser having a square-lattice structure. In this study we perform the 3D-FDTD calculation for the structure of an actual fabricated device. The device is based on band-edge resonance, and four band edges are present at the corresponding band edge point. For these band edges, we calculate the quality (*Q*) factor. The results show that the *Q* factor of a resonant mode labeled A_{1} is larger than that of other resonant modes; that is, lasing occurs easily in mode A_{1}. The device can thus achieve single-mode lasing oscillation. To increase the *Q* factor, we also consider the optimization of device parameters. The results provide important guidelines for device fabrication.

© 2005 Optical Society of America

## 1. Introduction

Recently, many novel applications using the photonic crystal (PC) have been proposed[1–6]. Among these applications, a unique semiconductor laser that uses band edge engineering to create photonic band structure has attracted a lot of attention[7–16]. In such a device, standing waves formed at each band edge are used as laser cavities. In previous studies, we reported a surface-emitting laser having a two-dimensional (2D) PC, by utilizing band-edges. Single longitudinal and/or lateral mode oscillation was successfully achieved over a large area[7]. We also demonstrated polarization mode control in a 2D square-lattice PC laser by controlling the geometry of the unit cell structure[9,10]. In these studies, we utilize the folded Γ point toward the Γ-X direction (Γ2 point). This band edge induces not only in-plane optical coupling, but also diffraction in the direction normal to the PC plane. As a result, the device functions as a surface-emitting laser. In the case of the square-lattice PC having a circular unit cell structure discussed in this paper, four band edges are present at the band edge of the Γ2 point. In ref [11], we reported the differences in features such as the mode profile of resonant mode, far field patterns, and quality (*Q*) factors that occur among the resonant modes originating in these band edge resonances for the case of a 2D square lattice PC slab, formed from a dielectric slab having circular air rods and with air cladding layers above and below. We discerned obvious differences in the mode profiles, far field patterns, and *Q* factors among these resonant modes. However, the actual device that we fabricated is formed by an active layer sandwiched by p- and n-cladding layers. The 2D PC structure is embedded in the cladding layer near the active layer.

In this study, we report characteristics of resonant modes, focusing mainly on the *Q* factor at each band edge of the Γ2 point for a structure similar to that of the actual device. Threshold gain is proportional to *Q* factor, and high *Q* factor is advantageous for achieving low-threshold lasing action. We use the 3D finite-difference time-domain (FDTD) method to calculate the mode profile and *Q* factor. In section 2, we describe the method for calculating the resonant mode. In section 3, we describe in detail the results of resonant mode calculations. We prove that *Q* factors differ depending on the band edges of the Γ2 point and especially on the unit cell structure size or air filling factor. Therefore, there exists an optimum unit cell structure size for achieving a low threshold current.

## 2. Calculation model and method

#### 2.1 Calculation of realistic structure

In this section, we explain the calculation method. We use the FDTD method [17] to calculate the mode profile and *Q* factor. The FDTD calculation method is almost the same as that employed in our previous study of a PC slab[11]. However, unlike the PC slab structure discussed in the previous study, the structure discussed in this study leads to some problems in calculating the mode profile and *Q* factor. One problem is that the effective refractive index difference of the PC of the structure discussed in this study is very small as compared with that of the PC slab. This is because the PC structure exists at the cladding layer, and only the evanescent component, which is leaked from the active layer, is influenced by modulation of the refractive index. To achieve a high *Q* factor, we must sufficiently increase the period of the PC. Furthermore, for the structure discussed in this study, light confinement at the active layer is very weak as compared with the case of the PC slab, because the refractive index difference between the active layer and the cladding layer is very small. Therefore, in order to attenuate the evanescent component sufficiently, the cladding layer must be of sufficient thickness. As a result, FDTD calculation must be performed for a very large model. However, obtaining the solution for such a large model requires very high computation power, very large memory size, and a very long time. This is the first problem. The second problem is as follows. We consider the 2D square-lattice PC, which has a circular unit cell structure. As mentioned previously, in the case of the square-lattice PC, four band edges are present at the band edge of the Γ2 point. For the structure discussed in this study, resonant frequency differences among these resonant modes are very small, because, as described previously, the effective refractive index difference is very small. Therefore, multiple resonant modes are excited simultaneously, making it difficult to excite only a specific resonant mode. Also, as mentioned later, one of the resonant modes has a *Q* factor much smaller than those of the other resonant modes. The resonant mode with a small *Q* factor is especially difficult to excite. This is the second problem. In order to avoid these two difficulties, we utilize symmetry of structure and resonant modes[15,16,18].

#### 2.2 Symmetry analysis

Figure 1 depicts the magnetic field distribution normal to the PC plane, obtained by the 2D plane wave (PW) expansion method using the effective refractive index[10,19]. The method for determining the effective index follows Ref.[8]. Thick black circles indicate the locations and shapes of lattice points. We can classify the resonant modes originating from the four band edges at the Γ2 point into A1, B1, and E modes. This nomenclature is based on the representation under group theory, which indicates the symmetry of the eigen mode at each band edge when the period of the PC is infinite. The characteristics of each mode are as follows. When we take the *x* and *y* axes as shown in Fig. 2, the structure is invariant against a symmetry operation in which *x* is changed to -*x*; i.e., is symmetrical with respect to the *y* axis. This symmetry operation is called σ_{x}. Similarly, the symmetry operations in which *y* is changed to -*y*, (*x*, *y*) to (*y*, *x*), and (*x*, *y*) to (-*y*, -*x*) are called σ
_{y}
, σ’
_{d}
, σ”
_{d}
, respectively. The structure is also invariant against rotations by 90, 180, and 270 degrees; such symmetry operations are called C_{4}, ${{\mathrm{C}}_{4}}^{2}$, and ${{\mathrm{C}}_{4}}^{-1}$, respectively. Thus, the square lattice structure with a circular unit cell structure is invariant against symmetry operations {E, C_{4}, ${{\mathrm{C}}_{4}}^{-1}$, ${{\mathrm{C}}_{4}}^{2}$, σ
_{x}
, σ
_{y}
, σ’
_{d}
, σ”
_{d}
}, where E is the symmetry operation that maintains the structure as it is. The magnetic field (H_{z}) distribution of the A_{1} mode is invariant against all symmetry operations, including E, C_{4}, ${{\mathrm{C}}_{4}}^{-1}$, and the others. On the other hand, the H_{z} distribution of the B_{1} mode is invariant against symmetry operations E, ${{\mathrm{C}}_{4}}^{2}$, σ
_{x}
, and σ
_{y}
, but not against symmetry operations C_{4}, ${{\mathrm{C}}_{4}}^{-1}$, σ’
_{d}
, σ”
_{d}
, which change the sign of Hz. Moreover, even the E mode is relatively complicated. The important point is that one of the E modes is a replica of the other given by 90 degrees rotation; these two modes are essentially the same and cannot be distinguished from each other, and thus are doubly degenerated[20,21].

#### 2.3 FDTD calculation method

Figure 3 shows the calculation model. As shown in Fig. 3, we divide the *x*-*y* plane into four parts by symmetrical planes; i.e., the two symmetrical planes A and B. Then we introduce the PMC (Perfect Magnetic Conductor) or PEC (Perfect Electric Conductor) boundary condition to these two planes. By setting the PEC boundary condition for both the A and B planes, we can excite the A_{1} mode or the B_{1} mode. By setting the PEC and PMC boundary conditions for the A and B planes, respectively, we can excite one of the doubly degenerated E modes. Conversely, by setting the PMC and PEC boundary conditions for the A and B planes, respectively, we can excite the other E mode. However, these two E modes are equivalent and exhibit no differences. Therefore, it is sufficient to investigate either of these E modes. In order to distinguish the A_{1} and B_{1} modes, we must introduce a symmetrical plane at an angle of 45 degrees to the symmetrical planes A and B. However, introducing such a symmetrical plane in FDTD calculation is difficult. Accordingly, we distinguish the A_{1} and B_{1} modes by changing the excitation points and excitation frequency between the A_{1} and B_{1} modes[11]. As mentioned later, because the *Q* factor is almost identical for the A_{1} and B_{1} modes and is sufficiently high, we can use this method to distinguish these two resonant modes. However, in order to investigate the characteristics of the E mode, which has a *Q* factor drastically lower than those of the A_{1} and B_{1} modes, it is convenient to completely differentiate between the A_{1} and B_{1} modes and the E mode. The 2D-PC structure of the experimental device is embedded on either side; i.e., in the n- and p-cladding layers. Such an asymmetrical structure with respect to the *z* direction generates an asymmetry mode along the *z* direction and makes analysis difficult. In order to avoid this, we render the calculation model a symmetrical structure with respect to the *z* direction and introduce the symmetrical boundary condition at the center of the active layer. As a result, the 2D-PC structure exists at both sides of the cladding layer. By adopting this model, we can distinguish a TE-like mode and a TM-like mode. Because the laser oscillation of the experimental device occurs in a TE-like mode, in this study we discuss only the TE-like mode. By introducing these symmetrical boundary conditions, we can reduce the calculation burden to one-eighth and can greatly reduce calculation time. In this study, we use these symmetry characteristics in calculating the resonant frequency, mode profile, and *Q* factor by 3D-FDTD simulation.

Here we introduce two important parameters concerning PC size. One is the number of air rods in the Γ-X direction (*N*), as shown in Fig. 3. Because we introduce symmetrical boundary conditions, the total number of air rods is equal to 2*N*×2*N*. The other is the air filling factor (*F*), which is defined as the area fraction occupied by the air rod per unit cell. The parameters used in the calculation model are a dielectric constant of 11.5, and an active layer thickness of 0.24 µm (or 0.6*a*, where *a* is the lattice constant and is equal to 0.4 µm.). The cladding layer has a dielectric constant of 10.3 and a thickness of 1.44 µm. The thickness of the air layer (*ε*=1) outside the cladding layer is 1.08 µm. The PC layer is placed 0.04 µm from the active layer, and has a thickness of 0.4 µm. The PC layer contains circular air rods (*ε*=1). The air rods are arranged in a square region as shown in Fig. 3. Mur’s second-order absorbing boundary condition is employed[22]. To estimate the resonant frequency at the corresponding band edge, the photonic band structure is also calculated by the 3D-FDTD method. For this calculation, Bloch’s and Mur’s second-order absorbing boundary conditions are employed for the horizontal and vertical boundary conditions, respectively. The parameters used in the FDTD calculation are Δ
_{x}
=Δ
_{y}
=Δ
_{z}
=1/10*a*, and Δ*t*=0.5Δ*x*/*c*, where *c* is the speed of light in a vacuum. The excitation method and the calculation method of the *Q* factor described in Ref[11] are employed. To separate the loss of the guided mode in the PC plane and the out-of-plane radiation loss in the direction normal to the PC plane, we decompose the *Q* factor (*Q*
_{total}) into a horizontal *Q* factor (*Q _{‖}*) and a vertical

*Q*factor (

*Q*

_{⊥})[23].

## 3. Calculation results

This section describes the results of calculation, including the resonant frequency, mode profile, and *Q* factors of the resonant modes. First, we show the resonant field distributions in Fig. 4. Figure 4(a) shows the magnetic field distributions normal to the PC plane at the center of the active layer in the A_{1} mode. The number of air rods in the Γ-X direction (*N*) is 50, and the air filling factor (*F*) is 10%. The rectangular region shown by the black dotted line in Fig. 4(a) indicates the PC region. Figure 4(a) shows that the magnetic field at the active layer is subject to strong feedback by the PC embedded at the cladding layer. The top-right inset of Fig. 4(a) shows a magnification of the field distribution near the center of the PC region. The field distribution indicated by the inset is consistent with that calculated by the 2D PW expansion method as shown in Fig. 1(a), indicating that the resonant field distributions are caused by band edge resonance. Figure 4(b) is obtained by the Fourier transform of the magnetic field of the PC area in Fig. 4(a). We can consider Fig. 4(b) as representing the *k*(wave number)-space pattern of the A_{1} mode. The white dotted square in Fig. 4(b) represents the first Brillouin zone. The *k*-space pattern of the A_{1} mode, as shown in Fig. 4(b), indicates that this mode is caused by the resonance at the Γ point. This is caused by finiteness of the structure and influence of the window function; the wave number is not fixed securely at the Γ point and spreads over a certain region around the Γ point. These results are in agreement with those obtained in a previous study using the PC slab[11]. We also confirm that other resonant modes such as the B_{1} and E modes, which are not shown here, are consistent with results obtained by study of PC slab and/or the 2D PW expansion method.

Next, in Fig. 5 we show the resonant frequency of each resonant mode as a function of *F* when *N*=50. The maximum value of *F* is 70% in Fig. 5, because the filling factor determined by the area of a circle inscribed in a square is 79.6%. Figure 5 shows the height of the resonant among the A_{1}, B_{1}, and E modes changes as a function of *F*. We can classify air filling factor into three regions depending on the relative height of the resonant frequency. When *F* is less than about 30%, the resonant frequency increases in the order of A_{1}, B_{1}, and E. We call this region A. When *F* is around 30%, the resonant frequency increases in the order of A_{1}, E, and B_{1}. We call this region B. In this region, the differences in resonant frequency among the resonant modes are very small. When *F* is more than about 30%, the resonant frequency increases in the order of E, A_{1}, and B_{1}. We call this region C.

We can obtain the same results by the 2D PW expansion method using the effective refractive index. In Fig. 6, we show the typical band structure of each region around the Γ2 point, as calculated by the 2D PW expansion method. From Fig. 6, we can see changes in not only the relative height of the resonant frequency, but also the band shape. For example, the shape of the band around the band edge A_{1} is upwards convex when *F* is less then 30%; that is, in region A or B, but becomes flatter toward Γ-X when *F* is more than 30%; that is, region C. Meanwhile, the shape of band around the band edge B_{1} is flat toward the Γ-X direction in region A, but is downwards convex in regions B and C. An interesting observation is that the boundary points between region A and region B and between region B and region C are triply degenerated.

Next, we explain the relationship between the stability of the guided mode and *F*. A large air filling factor is equivalent to a low equivalent refractive index for the PC layer. In the example discussed in this study, the PC layer exists at the cladding layer on both sides of the active layer and is located very close to the active layer. When the refractive index of the PC layer becomes low, light cannot be stably confined to the active layer and light leaks into the cladding layer. However, because an air layer exists outside the cladding layer, the guided mode is maintained by total reflection between the air and cladding layers and the mode cannot operate as an anti-guide. We explain this phenomenon by reference to Fig. 7, which shows the *k*-space patterns of the A_{1} mode when *F* is 10% or 60%. Considering the propagation of light normal to the PC plane, these *k*-space patterns are obtained by Fourier transformation of the electric field (E
_{x}
) at the center of the active layer, unlike the case in Fig. 4(b). In Fig. 7, the circles indicated by a solid line, a broken line, and a dotted line are the light cones determined by the refractive index of the active layer, the cladding layer, and air, respectively. Figure 7(a) reveals that when *F* is 10%; that is, when the resonant frequency is 0.304, only the wave number components around *k*=(0, 0) exist inside the air light cone. These wave number components are the origin of the surface-emitting component. The wave number components around *k*=(0, 2*π*/*a*) and/or *k*=(0,-2*π*/*a*) exist between the active layer light cone and the cladding layer light cone. This means that these wave number components are confined to the active layer by total reflection between the active and cladding layers and that they generate a guided mode. On the other hand, Fig. 7(b) shows that when *F* is 60%; that is, when the resonant frequency is 0.318, the wave number components around *k*=(0, 2*π*/*a*) and/or *k*=(0, -2*π*/*a*) exist inside the cladding light cone but outside the air light cone. Thus, these wave number components are confined by total reflection between the cladding and air layers. Even in this case, the wave number components existing inside the air light cone are components around *k*=(0, 0). Therefore, the wave number components around *k*=(0, 2*π*/*a*) and/or *k*=(0, -2*π*/*a*) cannot be emitted to air. However, this causes a decrease in light confinement factor at the active layer and an accompanying drop in *Q* factor. As a result, the threshold current is expected to increase. The horizontal line at the frequency of 0.3116 in Fig. 5 indicates the threshold above which light leaks to the cladding layer. Thus, the air filling factor greatly influences the stability of the guided mode and the light confinement factor. This phenomenon is caused by the 2D-PC layer at both sides of the cladding layer, with the effective refractive index becoming low. When the PC layer is embedded on either side; i.e., in the n- and p-cladding layers as in our experimental device or in a PC slab composed of three simple layers, this phenomenon does not occur.

Figure 8 shows the relationship between the total quality factor (*Q*
_{total}) and *F* for each resonant mode at *N*=50. In our previous study, we reported that *Q*
_{total} of the A_{1} mode achieves over 5000 at a PC size of 50×50, but in the structure discussed in this study *Q*
_{total} of the A_{1} mode is about 1400 at a PC size of 100×100 (*N*=50). Thus, the effect of the PC structure can be considered smaller for this structure than for the PC slab structure. Consequently, to achieve a sufficiently high *Q* factor to lase, we need a larger PC. In fact, our experimental device has a PC size of 1000×1000. The structure discussed in this study has a PC size of 100×100 at the largest and is considerably smaller than the experimental device, but we still consider the model to be useful for analyzing the *Q* factor tendency with different resonant modes and filling factors. From Fig. 8, the height of *Q*
_{total} between the A_{1} and B_{1} modes is inverted when *F* is over 50%, but when *F* is less than about 50%, *Q*
_{total} decreases in the order of A_{1}, B_{1}, and E. This result shows that among the corresponding resonant modes, the A_{1} mode allows lasing to occur relatively easily. Meanwhile, *Q*
_{total} of the E mode is considerably smaller than those of the other resonant modes. These results are in agreement with a previous study of a PC slab. We can obtain the maximum *Q*
_{toal} at *F*=10%. To achieve a low threshold current, a device should be fabricated with this filling factor. Recently, we experimentally confirmed that the threshold current decreases at an air filling factor of around 10%[24].

To consider the *Q* factor characteristics as a function of *F*, we decompose *Q*
_{total} into *Q _{‖}* and

*Q*

_{⊥}, as shown in Fig. 9, which shows

*Q*and

_{‖}*Q*

_{⊥}for each resonant mode plotted against

*F*. Figure 9 shows that, in the A

_{1}and B

_{1}modes,

*Q*

_{total}is almost entirely determined by

*Q*;

_{‖}*Q*is much smaller than

_{‖}*Q*

_{⊥}. On the other hand, in the E mode,

*Q*

_{total}is almost entirely determined by

*Q*

_{⊥};

*Q*

_{⊥}is much more smaller than

*Q*. The reason is that in a PC of infinite size, the A

_{‖}_{1}and B

_{1}modes do not couple to free space normal to the PC plane, because of a symmetry mismatch, but the E mode does couple to free space normal to the PC plane[25,26]. The trend of

*Q*with

_{‖}*F*is notable. For all modes,

*Q*exhibits local maximums at around

_{‖}*F*=10% and

*F*=45%, and a minimum at around

*F*=30%. As a result, the graph of

*Q*against

_{‖}*F*shows dual peaks. This characteristic is similar to that of the coupling coefficient

*κ*of a 2nd order 1D DFB laser, which shows dual peaks with a node at

*F*=50%[27]. It is also similar to the results derived from Fig. 5 or Fig. 6.

Figure 10 shows the frequency difference between the maximum and minimum resonant frequencies; in other words, the bandgap width at the Γ2 point, with changing *F*. The solid line is obtained by the 3D-FDTD calculation shown in Fig. 5 and the dotted line is obtained by the 2D PW expansion method using the effective refractive index shown in Fig. 6. The absolute values are not consistent between the 3D-FDTD and 2D PW expansions, but the tendencies are consistent. An important point is that the trends observed in Fig. 10 are similar to those for *Q _{‖}*. If we consider this frequency difference to be the stopband width, light receives a large feedback when the frequency difference is large; that is,

*Q*becomes large. This result indicates that we can easily estimate the

_{‖}*Q*characteristics from the resonant frequency difference. As shown in Fig. 6, we can estimate the resonant frequency difference from the 2D PW expansion method using the effective refractive index. The 3D-FDTD calculation, which requires a very long time and very high computation power, is not required. Even more convenient is that the 2D-FDTD calculation is also unnecessary, and we can estimate the

_{‖}*Q*characteristics simply by using the 2D PW expansion method. As mentioned before,

_{‖}*Q*

_{total}is almost entirely determined by

*Q*for the A

_{‖}_{1}mode; therefore, we can regard the characteristics of

*Q*to be those of

_{‖}*Q*

_{total}. Because

*Q*

_{total}of the A

_{1}mode is larger than that of other modes, the A

_{1}mode allows relatively easy lasing. As a result, we can obtain the optimum filling factor of the device by this estimation. However, it is also important to consider application of this estimation method. For example, in the case of a PC slab which has a strong refractive index modulation, estimating the effective refractive index and effective index difference is very difficult, leading to difficulty in determining the resonant frequency by the 2D PW expansion method. In summary, a significant point to consider is whether it is reasonable to perform the 2D PW expansion method using the effective refractive index. We consider this estimation to be valid, as long as it is possible to replace the 3D model with the 2D model. This may prove very difficult and the results of future investigations in this area will be reported elsewhere.

## 4. Summary

The optical properties of a two-dimensional (2D) square-lattice photonic crystal (PC) surface-emitting laser have been investigated in detail using the 3D-FDTD method. The structure studied by calculation was selected to be similar to that of our experimental device composed of an n-cladding layer, a p-cladding layer, and an active layer. The PC is embedded at the cladding layer and is close to the active layer. First, we show the field distributions of the resonant mode, and confirm that the obtained field distribution is caused by resonance at the Γ point. Then we investigate the resonant frequency of each resonant mode with changing air filling factor. The stability of the guided mode is confirmed to be related to the air filling factor with the possibility of instability in the guided mode arising when the air filling factor becomes large. Thus, the air filling factor greatly influences the stability of the guided mode and the light confinement factor. This phenomenon arises from the 2D-PC layer existing on both sides of the cladding layer and having a low effective refractive index. This phenomenon does not occur when the PC layer is embedded on either side of the n- and p-cladding layers as in our experimental device or in a PC slab which is composed of three simple layers.

Next, we investigate the *Q* factors of the resonant modes. *Q* calculations indicate that *Q*
_{total} of the A_{1} mode is about 1400 at a PC size of 100×100 (*N*=50). Compared with the result for the PC slab, *Q*
_{total} obtained by this study is considerably smaller, indicating that, in order to achieve a sufficiently high *Q* factor to lase, our PC must be larger than the PC slab. Among the resonant modes, the A_{1} mode has the highest *Q*
_{total}; that is, among the resonant modes, the A_{1} mode allows relatively easy lasing. The A_{1} mode is shown to attain a maximum *Q*
_{total} when the air filling factor is 10%. Results of our recent experiments indicate that the threshold current decreases at an air filling factor of around 10%. The characteristics of *Q*
_{total} are clarified by decomposing *Q*
_{total} into *Q _{‖}* and

*Q*

_{⊥}. This decomposition shows that

*Q*

_{total}is almost entirely determined by

*Q*for the A

_{‖}_{1}and B

_{1}modes, but that of the E mode is almost entirely determined by

*Q*

_{⊥}. The A

_{1}and B1 modes do not easily couple to free space normal to the PC plane because of the symmetry mismatch, but the E mode can couple to free space normal to the PC plane.

Finally, the existence of dual peaks in the graph of *Q _{‖}* against F is considered for all resonant modes. This result is similar to that for the resonant frequency difference between the maximum and minimum values of frequency at each resonant mode or bandgap width at the Γ2 point. Thus, the band gap width at the Γ2 point can be regarded to be the stopband width. The influence of bandgap width at the Γ2 point on

*Q*factor may be estimated by the 2D PW expansion method, rather than using the complex 3D-FDTD calculation, allowing the estimation to be performed simply and conveniently.

## References and links

**1. **E. Yablonovitch, “Inhibited Spontaneous Emission in Solid-State Physics and Electronics,” Phys. Rev. Lett. **58**, 2059–2062 (1987). [CrossRef] [PubMed]

**2. **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]

**3. **S. Noda, K. Tomoda, N. Yamamoto, and A. Chutinan, “Full Three-Dimensional Photonic Crystals at Near-Infrared Wavelengths,” Science **289**, 604–606 (2000). [CrossRef] [PubMed]

**4. **S. Noda, M. Imada, and A. Chutinan, “Trapping and emission of photons by a single defect in a photonic bandgap structure,” Nature **407**, 608–610 (2000). [CrossRef] [PubMed]

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

**6. **S. Noda and T. Baba, Eds., *Roadmap on Photonic Crystals*, (Kluwer Academic, New York, 2003).

**7. **M. Imada, S. Noda, A. Chutinan, T. Tokuda, M. Murata, and G. Sasaki, “Coherent two-dimensional lasing action in surface-emitting laser with triangular-lattice photonic crystal structure,” Appl. Phys. Lett. **75**, 316–318 (1999). [CrossRef]

**8. **M. Imada, A. Chutinan, S. Noda, and M. Mochizuki, “Multidirectionally distributed feedback photonic crystal lasers,” Phys. Rev. B **65**, 195306 (2002). [CrossRef]

**9. **S. Noda, M. Yokoyama, M. Imada, A. Chutinan, and M. Mochizuki, “Polarization mode control of two-dimensional photonic crystal laser by unit cell structure design,” Science **293**, 1123–1125 (2001). [CrossRef] [PubMed]

**10. **M. Yokoyama and S. Noda, “Polarization mode control of two-dimensional photonic crystal laser having a square lattice structure,” IEEE J. Quantum Electron. **39**, 1074–1080 (2003). [CrossRef]

**11. **M. Yokoyama and S. Noda, “Finite-Difference Time-Domain Simulation of Two-Dimensional Photonic Crystal Surface-Emitting Laser having a Square-Lattice Slab Structure,” IEICE Trans. Electron. **E87-C**, 386–392 (2004).

**12. **M. Meier, A. Mekis, A. Dodabalapur, A. Timko, R. E. Slusher, J. D. Joannopoulos, and O. Nalamasu, “Laser action from two-dimensional distributed feedback in photonic crystals,” Appl. Phys. Lett. **74**, 7–9 (1999). [CrossRef]

**13. **K. Inoue, M. Sasada, J. Kawamata, K. Sakoda, and J. W. Haus, “A Two-Dimensional Photonic Crystal Laser,” Jpn. J. Appl. Phys. **38**, L157–L159 (1999). [CrossRef]

**14. **M. Meier, A. Dodabalapur, J. A. Rogers, R. E. Slusher, A. Mekis, A. Timko, C. A. Murray, R. Ruel, and O. Nalamasu, “Emission characteristics of two-dimensional organic photonic crystal lasers fabricated by replica molding,” J. Appl. Phys. **86**, 3502–3507 (1999). [CrossRef]

**15. **R. Colombelli, K. Srinivasan, M. Troccoli, O. Painter, C. Gmachl, D. M. Tennant, A. M. Sergent, D. L. Sivco, A. Y. Cho, and F. Capasso, “Quantum cascade surface-emitting photonic crystal laser,” Science **302**, 1374–1377 (2003). [CrossRef] [PubMed]

**16. **K. Srinivasan, O. Painter, R. Colombelli, C. Gmachl, D.M. Tennant, A.M. Sergent, D.L. Sivco, A.Y. Cho, M. Troccoli, and F. Capasso, “Lasing mode pattern of a quantum cascade photonic crystal surface-emitting microcavity laser,” App. Phys. Lett. **84**, 4164–4166 (2004). [CrossRef]

**17. **K. S. Yee, “Numerical Solution of Initial Boundary Value Problem Involving Maxwell’s Equations in Isotropic Media,” in *Proceedings of IEEE Conference on Antennas and Propagat. AP-14* (Institute of Electrical and Electronics Engineers, New York, 1966), pp. 302–307.

**18. **M. Okano and S. Noda, “Analysis of multimode point-defect cavities in three-dimensional photonic crystals using group theory in frequency and time domains,” Phys. Rev. B **70**, 125105 (2004). [CrossRef]

**19. **M. Plihal, A. Shambrook, and A. A. Maradudin, “Two-dimensional photonic band structures,” Opt. Comm. **80**, 199–204 (1991). [CrossRef]

**20. **K. Sakoda, “Symmetry, degeneracy, and uncoupled modes in two-dimensional photonic lattices,” Phys. Rev. B **52**, 7982–7986 (1995). [CrossRef]

**21. **K. Sakoda, *Optical Properties of Photonic Crystals*, (Springer Verlag, Berlin, 2001).

**22. **G. Mur, “Absorbing Boundary Conditions for the Finite-Difference Approximation of the Time-Domain Electromagnetic-Field Equations,” in *Proceedings of IEEE Conference on Electromagn. Compat. EMC-23* (Institute of Electrical and Electronics Engineers, New York, 1981), pp. 377–382.

**23. **O. J. Painter, J. Vuckovic, and A. Scherer, “Defect modes of a two-dimensional photonic crystal in an optically thin dielectric slab,” J. Opt. Soc. Am. B **16**, 275–285 (1999). [CrossRef]

**24. **D. Ohnishi, T. Okano, M. Imada, and S. Noda, “Room temperature continuous wave operation of a surface-emitting two-dimensional photonic crystal diode laser,” Opt. Express **12**, 1562–1568 (2004). [CrossRef] [PubMed]

**25. **T. Ochiai and K. Sakoda, “Dispersion relation and optical transmittance of a hexagonal photonic crystal slab,” Phys. Rev. B **63**, 125107 (2001). [CrossRef]

**26. **S. Fan and J. D. Joannopoulos, “Analysis of guided resonances in photonic crystal slabs,” Phys. Rev. B **65**, 235112 (2002). [CrossRef]

**27. **M. Imada, S. Noda, H. Kobayashi, and G. Sasaki, “Characterization of a Distributed Feedback Laser with Air/Semiconductor Gratings Embedded by the Wafer Fusion Technique,” IEEE J. Quantum Electron. **35**, 1277–1283 (1999). [CrossRef]