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₄, ${{\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₄, ${{\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₁ mode is invariant against all symmetry operations, including E, C₄, ${{\mathrm{C}}_4}^{-1}$, and the others. On the other hand, the H_z distribution of the B₁ mode is invariant against symmetry operations E, ${{\mathrm{C}}_4}^2$, σ _x , and σ _y , but not against symmetry operations C₄, ${{\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].

 

Fig. 1. (a), (b), (c), (d) Magnetic field distributions at band edges A₁, B₁, E, and E, respectively. The amplitudes of the magnetic fields in the direction perpendicular to the plane are indicated by white and black areas, which denote positive and negative regions, respectively. The thick black circles indicate the shapes and locations of lattice points.

Download Full Size | PDF

 

Fig. 2. Coordinates and symmetry of the square lattice structure having a circular unit cell structure

Download Full Size | PDF

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₁ mode or the B₁ 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₁ and B₁ 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₁ and B₁ modes by changing the excitation points and excitation frequency between the A₁ and B₁ modes[11]. As mentioned later, because the Q factor is almost identical for the A₁ and B₁ 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₁ and B₁ modes, it is convenient to completely differentiate between the A₁ and B₁ 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 2N×2N. 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.6a, 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/10a, 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].

 

Fig. 3. Schematic of 3D-FDTD calculations. The device is composed of an n-cladding layer, a p-cladding layer, and an active layer. The 2D PC is embedded at the cladding layer close to the active layer.

Download Full Size | PDF

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₁ 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₁ mode. The white dotted square in Fig. 4(b) represents the first Brillouin zone. The k-space pattern of the A₁ 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₁ 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.

 

Fig. 4. (a) Magnetic field distributions of A₁ mode normal to the photonic crystal plane at the center of the active layer. The PC size is 100×100. (b) Wave number space (k-space) patterns of (a). The figure is obtained by Fourier transformation of the magnetic field of the photonic crystal area in Fig 4(a). The white dotted square represents the first Brillouin zone.

Download Full Size | PDF

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₁, B₁, 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₁, B₁, and E. We call this region A. When F is around 30%, the resonant frequency increases in the order of A₁, E, and B₁. 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₁, and B₁. We call this region C.

 

Fig. 5. The resonant frequency of each resonant mode versus air filling factor. We can classify air filling factor into three regions depending on the by the relative height the resonant frequency. The horizontal line at a frequency of 0.3116 is the cladding light line, above which light leaks to the cladding layer.

Download Full Size | PDF

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₁ 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₁ 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.

 

Fig. 6. (a) (b) (c) Typical band structures of regions A, B, and C around the Γ2 point, respectively. These band structures are calculated by the 2D PW expansion method. The dielectric constants inside the unit cell are 9.26, 9.86, and 10.01, respectively. The dielectric constants outside the unit cell are 10.60, 10.47, and 10.38 respectively. The air filling factors are 8.50%, 30.0%, and 50.0%, respectively.

Download Full Size | PDF

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₁ 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.

 

Fig. 7. (a)(b) The wave number space (k-space) patterns of the electric field (E _x ) of the A₁ mode, when the air filling factors are 10% and 60%, respectively. 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.

Download Full Size | PDF

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₁ mode achieves over 5000 at a PC size of 50×50, but in the structure discussed in this study Q _total of the A₁ 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₁ and B₁ modes is inverted when F is over 50%, but when F is less than about 50%, Q _total decreases in the order of A₁, B₁, and E. This result shows that among the corresponding resonant modes, the A₁ 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₁ and B₁ 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₁ and B₁ 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.

 

Fig. 8. Q factors as functions of air filling factor for A₁, B₁, and E modes. Size of the PC is 100×100. The Q factor of the A₁ mode is larger than those of the other resonant modes. This result indicates that, among the resonant modes, the A₁ mode allows relatively easy lasing.

Download Full Size | PDF

 

Fig. 9. (a)(b) Horizontal Q factor (Q_‖) and vertical Q factor (Q _⊥) as functions of the air filling factor for A₁, B₁, and E modes. Q_‖ and Q _⊥ represent losses of guided mode and out-of-plane radiation, respectively, in the direction normal to the photonic crystal plane.

Download Full Size | PDF

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₁ mode; therefore, we can regard the characteristics of Q_‖ to be those of Q _total. Because Q _total of the A₁ mode is larger than that of other modes, the A₁ 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.

 

Fig. 10. Frequency differences between maximum and minimum resonant frequencies for each of the A₁, B₁, and E resonant modes. The solid line is obtained by the 3D-FDTD calculation and the dotted line is obtained by the 2D PW expansion method using the effective refractive index. The characteristic trend is similar to that for Q_‖.

Download Full Size | PDF

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₁ 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₁ mode has the highest Q _total; that is, among the resonant modes, the A₁ mode allows relatively easy lasing. The A₁ 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₁ and B₁ modes, but that of the E mode is almost entirely determined by Q _⊥. The A₁ 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.