## Abstract

Three-dimensional coupled-wave theory is extended to model triangular-lattice photonic-crystal surface-emitting lasers with transverse-electric polarization. A generalized coupled-wave equation is derived to describe the sixfold symmetry of the eigenmodes in a triangular lattice. The extended theory includes the effects of both surface radiation and in-plane losses in a finite-size laser structure. Modal properties of interest including the band structure, radiation constant, threshold gain, field intensity profile, and far-field pattern (FFP) are calculated. The calculated band structure and FFP, as well as the predicted lasing mode, agree well with experimental observations. The effect of air-hole size on mode selection is also studied and confirmed by experiment.

© 2013 Optical Society of America

## 1. Introduction

Two-dimensional (2D) photonic-crystal surface-emitting lasers (PC-SELs) are becoming increasingly important owing to their promising functionality and improved performance compared to conventional semiconductor lasers. A number of successful demonstrations underpinning this promise have already been conducted [1–13]. By utilizing the band edge of the photonic band structure, single longitudinal and transverse mode oscillation in two dimensions has been achieved with a large lasing area, which enables high-power, single-mode operation [5,6]. The output beam of such devices is generally emitted in the direction normal to the 2D PC plane and has a small beam divergence angle of less than 1° owing to the large area of coherent oscillation [5, 6]. Furthermore, both the polarization and beam pattern [2, 6, 12, 13], as well as the beam direction [11], of the laser output can be controlled by appropriate design of the PC geometry. In addition, the lasing wavelengths of PC-SELs have spanned a wide range from near-infrared to blue-violet, mid-infrared, and terahertz regimes [6–10].

To further improve the device performance or exploit new functionalities of PC-SELs, it is important to develop an accurate first-principles model to describe the lasing properties in these structures. Realistic devices of PC-SELs require three-dimensional (3D) analysis, because both the surface emission in the vertical direction and the 2D optical feedback in the PC plane need to be taken into account. In addition, due to two important features of the PC laser cavities – a large area that is typically above 200 periods for one side length and increasingly complicated air-hole shape designs [2, 11, 21, 22] – it is difficult to model this type of laser using brute-force computer simulation methods such as the finite-difference time-domain (FDTD) [14, 15] or finite-element [10] methods.

One of the principal analytical tools that has been widely used to study laser devices incorporating periodic structures is the one-dimensional (1D) coupled-wave theory (CWT) originally proposed by Kogelnik and Shank in 1972 [16]. Recently, with the development of PC-SELs, this original theory has been extended, albeit specifically for a 2D model [4, 17–20]. More recently, we proposed a 3D CWT model that affords an exact analytical treatment of the full 3D structure and achieves a very accurate and efficient analysis of the band-edge modes of square-lattice PC-SEL devices [21–23]. Moreover, this 3D CWT framework has been extended in a fairly general fashion [24]. Now, it is not only applicable to the general centered-rectangular lattice PC structures [13] but also capable of treating the non-Γ-point modes in the vicinity of the second-order Γ point. However, all of our previous works on the 3D CWT have thus far been restricted to lattice structures for which the eigenstate of light can be described by four dominant wavevectors (i.e., basic waves). As a result, they are intrinsically applicable to lattice structure with *C*_{2v} or *C*_{4v} symmetry [25] and are unable to treat PCs with higher symmetry, e.g., triangular-lattice PCs with *C*_{6v} symmetry [19,25]. In this study, we further extend the 3D CWT works [21–24] to develop a coupled-wave model for triangular-lattice PC-SELs. We investigate various modal properties of interest, including the band structure, radiation constant, threshold gain, field intensity envelope profile, and far-field pattern (FFP). The calculated results of the band structure and FFP of the predicted lasing mode are compared with experimental data to validate our theory. Finally, the effect of the air-hole size on the mode selection of the laser device is studied and discussed.

The remainder of this paper is organized as follows. Section 2 describes the PC geometry and derivations of the coupled-wave equations. Section 3 presents numerical results, their comparison with experimental results, and discussions. Section 4 concludes this work by summarizing our findings.

## 2. Formulation

#### 2.1. Device geometry and photonic-crystal design

A schematic of the PC-SEL device investigated in this study is shown in Fig. 1(a). The PC layer is embedded inside a multilayer structure that is assumed to support only a single fundamental guided mode. Details of the device are described in [26], and the structural parameters are listed in Table 1. The average refractive index of the PC layer is given by
${n}_{av}=\sqrt{f{n}_{a}^{2}+\left(1-f\right){n}_{b}^{2}}$, where *n _{a}* is the refractive index of air,

*n*is the refractive index of the background dielectric material (GaAs), and

_{b}*f*is the air-hole filling factor (i.e., the fraction of the area of a unit cell occupied by air holes). For circular-shape air holes, $f=\left(2\pi /\sqrt{3}\right){r}^{2}/{a}^{2}$, where

*r*is the hole radius and

*a*is the lattice constant. A typical band structure calculated by the 2D plane-wave expansion method [3] is shown in Fig. 1(b). We are particularly interested in the band-edge modes in the vicinity of the second-order Γ point (red circle), where the group velocity of light becomes zero. At these band edges, the PC lattice serves both to provide the 2D distributed feedback and to couple the in-plane guided waves into the surface normal to construct the output.

Figure 2(a) shows a schematic of a triangular-lattice PC. The primitive translation vectors of the triangular lattice, **a _{1}** and

**a**, can be expressed as

_{2}*m*and

*n*are integers.

Light propagating inside a PC must obey Bloch’s theorem, and the Bloch wave state *u*(**r**) can be expressed as

*a*represents the field amplitude of a wave of the order of (

_{mn}*m*,

*n*) and Δ

**k**= Δ

_{x}β_{0}

*x̂*+ Δ

_{y}β_{0}

*ŷ*represents the wavenumber deviation from the Γ point [24]. Specifically, in the vicinity of the second-order Γ point (Δ

**k**≃ 0), only the field amplitudes of the waves with |

**G**

*| =*

_{mn}*β*

_{0}are dominant [3, 21, 24]. We refer to these waves as basic waves:

*R*

_{1},

*S*

_{1},

*R*

_{2},

*S*

_{2},

*R*

_{3}, and

*S*

_{3}, as illustrated in Fig. 2(b).

#### 2.2. Derivation of 3D coupled-wave equations

By following Refs. [21, 23, 24], we seek solutions to

**E**(

**r**) is

*e*,

^{iωt}*k*

_{0}(=

*ω/c*) is the free-space wavenumber,

*ω*is the angular frequency,

*c*is the velocity of light in free space, and

*ñ*is the refractive index (a complex number) satisfying ${k}_{0}^{2}{\tilde{n}}^{2}\left(\mathbf{r}\right)\simeq {k}_{0}^{2}{n}^{2}\left(\mathbf{r}\right)+2i{k}_{0}{n}_{0}\left(z\right)\tilde{\alpha}\left(z\right)$[23];

*n*

_{0}(

*z*) represents the average refractive index of the material at position

*z*, and

*α̃*(

*z*) represents the gain (

*α̃*> 0) or loss (

*α̃*< 0) in each region. For TE polarization, the electric field can be assumed as (

*E*,

_{x}*E*, 0) and expanded as

_{y}*m*≡

_{x,mn}*m*

_{0x}(

*m*,

*n*) + Δ

*and*

_{x}*n*=

_{y,mn}*n*

_{0y}(

*m*,

*n*) + Δ

*. Similarly, the periodic function of*

_{y}*n*

^{2}(

**r**) can be expanded as

*ξ*(

_{m,n}*z*) = 0 outside the PC layer, and air holes within the PC region are assumed to have perfectly vertical sidewalls [22]; thus, ${n}_{0}^{2}\left(z\right)$ and

*ξ*are constant within every region.

_{m,n}By substituting Eqs. (5) and (6) into Eq. (4) and collecting all terms that are multiplied by the factor *e*^{−imxβ0x−inyβ0y}, we obtain

By introducing a proper linear combination of *E _{x,m,n}*(

*z*) and

*E*(

_{y,m,n}*z*), i.e., (

*n*−

_{y}E_{x,m,n}*m*) as described in [21, 24], we obtain

_{x}E_{y,m,n}We are interested in the modes in the vicinity of the second-order Γ point, for which Δ* _{x}*, Δ

*≪ 1. We assume that all of the six basic waves have an identical vertical field profile Θ*

_{y}_{0}(

*z*) [24], which satisfies

*β*(

*β*≃

*β*

_{0}) and the vertical field profile Θ

_{0}(

*z*) can be calculated by employing the transfer matrix method (TMM) [21], and Θ

_{0}is normalized as ${\int}_{-\infty}^{\infty}{\mathrm{\Theta}}_{0}\left(z\right){\mathrm{\Theta}}_{0}^{*}\left(z\right)dz=1$. We further assume that basic waves are transverse waves. Thus, for basic waves, we can express the

*E*and

_{x,m,n}*E*as

_{y,m,n}*m,n*) ∈ {(1, 0), (−1, 0), (0, 1), (0, −1), (1, 1), (−1, −1)} which correspond, in order, to {

*R*

_{1},

*S*

_{1},

*R*

_{2},

*S*

_{2},

*R*

_{3},

*S*

_{3}} for

*A*(

_{m,n}*x*,

*y*), as shown in Fig. 2(b);

*ρ*and

_{mn}*η*denote the perturbation of the basic waves’ polarizations [24]. By substituting Eqs. (12) and (13) in Eq. (10), we obtain

_{mn}*δ*=

*β*−

*β*

_{0}=

*n*−

_{eff}ω/c*β*

_{0}is the deviation from the Bragg condition,

*n*is the effective refractive index for the fundamental guided mode,

_{eff}*α*is the modal loss defined by $\alpha =\frac{{k}_{0}}{{\beta}_{0}}{\int}_{-\infty}^{\infty}{n}_{0}\left(z\right)\tilde{\alpha}\left(z\right){\left|{\mathrm{\Theta}}_{0}\left(z\right)\right|}^{2}dz$. Note that the integrals in Eqs. (15) and (16) extend only over the PC region because

*ξ*= 0 outside that range. The field profiles [i.e.,

_{m,n}*E*(

_{x,m′,n′}*z*) and

*E*(

_{y,m′,n′}*z*)] of radiative (

*m′*=

*n′*= 0) and high-order (|

*m′*

^{2}+

*n′*

^{2}| > 1) waves can be expressed in terms of basic waves as described in Appendix A.

The 3D coupled-wave Eq. (15) derived above includes both the in-plane high-order coupling and vertical diffraction effects, neither of which was included or appropriately described in the previous 2D CWT models for triangular-lattice PC-SELs [4, 19, 27]. It is important to note that, Eq. (15) is a very general formula capable of treating any type of centered-rectangular lattice [13] by replacing the basic waves terms. Furthermore, the derivative terms of the basic waves, which are crucial factors for treating the finite-size laser structures but not considered in the previous 3D CWT model for centered-rectangular lattices [24], are now included.

## 3. Numerical results and discussion

In this section, we present numerical results of the coupled-wave equations derived above. By solving the coupled-wave Eq. (15), we can obtain various properties of interest, including the band structure, radiation constant, threshold gain, field intensity envelope profile within the device, and FFPs. To validate our CWT analysis, we also compare some of the calculated results with experimental observations. Unless otherwise noted in the following examples, the lattice constant *a* = 341 nm, the air-hole shape is circular with air-hole filling factor *f* = 0.15 (corresponding to an *r/a* ratio of 0.20, where *r* is the air-hole radius), and the structural parameters listed in Table 1 are used.

#### 3.1. Band structure and radiation constant

By considering an infinitely periodic structure, that is, by neglecting the derivative terms in the coupled-wave Eq. (15), the band structure in the vicinity of the second-order Γ point can be obtained. The calculated band structure in the two characteristic directions of the triangular lattice (Γ–*X* and Γ–*J* directions) is shown in Fig. 3(a), where the six eigenmodes are referred to as *A*, *B*_{1}, *B*_{2}, *C*, *D*_{1}, and *D*_{2}, respectively, in the order of increasing frequency. The leakiness of these eigenmodes is quantified in terms of the radiation constant *α _{r}* (

*α*= 2

_{r}*α*: the modal power loss due to the surface emission) [21, 24]. Figure 3(b) shows the calculated radiation constant plotted against the in-plane wavevector. It can be seen from the figure that modes

*D*

_{1}and

*D*

_{2}lead to a substantial radiation loss at the band edge (i.e., at the Γ point), whereas modes

*A*,

*B*

_{1},

*B*

_{2}, and

*C*are not leaky (

*α*= 0). This contrast is closely related with the symmetry of their field patterns with respect to the air holes. Figure 3(c) shows the calculated field patterns of the band-edge modes

*A*–

*D*

_{2}. We note that the electric fields of modes

*D*

_{1}and

*D*

_{2}are symmetric with respect to the air holes, and thus, they lead to significant loss. In contrast, the electric fields of modes

*A*,

*B*

_{1},

*B*

_{2}, and

*C*are antisymmetric, leading to zero radiation loss [3, 21]. As a consequence, the antisymmetric (nonradiative) modes are favored for lasing. The lasing mode is determined, as described in the following section, by taking into account the total losses (both the surface emission and in-plane losses) in a finite-size laser structure.

A noteworthy fact is that as the in-plane wavevector is slightly detuned from the band edge, most of the antisymmetric modes start to become significantly more lossy, regardless of the directions, thereby ensuring that band-edge modes keep lasing stably. However, the only exception is mode *C*, which has a slow-group-velocity region in the Γ–*J* directions. (We refer to modes in this region as flat-band modes.) Because the radiation constant becomes noticeably small within this region, the flat-band modes may also participate in lasing mode competition. This type of mode competition was reported in Ref. [28], where a multiple-mode lasing action occurring at both the band edge and the flat-band modes in the Γ–*J* direction was observed.

#### 3.2. Threshold gain and field intensity envelope

In this section, we solve the coupled-wave Eq. (15) in a finite-size system to evaluate the threshold gain of the eigenmodes. We discretize Eq. (15) by using the staggered-grid finite-difference method [23]. It is important to note here that, for the triangular-lattice PC configuration, it is crucial to both discretize the computational domain on a hexagonal grid [29] instead of the commonly used square grid and use appropriate boundary conditions (Appendix B). We focus our analysis on the antisymmetric band-edge modes mentioned above (i.e., modes *A*, *B*_{1}, *B*_{2}*,* and *C*), which are most likely to be lasing modes. In the following calculations, we assume an absorbing boundary condition and the radius of the circular-shape computational region *L* is set to be *L* = 30 *μ*m (Appendix B).

For finite-size structures, threshold gain and mode frequency correspond to the imaginary and real parts of the normalized eigenvalue, (*δ* + *iα*)*L*[23]. As an example, we show in Fig. 4(a) a plot of the normalized threshold gain (*αL*) as a function of deviation from the Bragg condition (*δL*). Though a large number of modes are found to exist as a result of the numerical calculation, we identify the fundamental modes (indicated by arrows) using the techniques described in Refs. [19, 23]. We note that there exists an additional mode that has a low threshold comparable to that of the fundamental band-edge modes. We refer to this mode as the *W* mode. The normalized threshold gains of the low-threshold modes indicated by arrows in Fig. 4(a) are listed in Table 2. As the lasing action occurs at the mode with the lowest threshold gain (loss), Table 2 indicates that mode *C* [indicated by red arrow in Fig. 4(a)] is favored for lasing with a large threshold-gain discrimination of over 20 cm^{−1} against the competing side modes (i.e., modes *B*_{1} and *B*_{2}). This fact reveals that stable single-mode operation at mode *C* can be possibly achieved for triangular-lattice PC-SELs.

To gain insight into the nature of the modes indicated in Fig. 4(a), it is useful to examine their field intensity envelopes throughout the entire laser cavity. The slowly-varying envelope profile (function) modulates the fast-varying Bloch waves in the PC lattice and can be determined by [20, 23]

*A*,

*B*

_{1},

*B*

_{2}, and

*C*) is relatively confined to the central part of the laser cavity, whereas mode

*W*exhibits a null intensity in the center. The exact physical interpretation of the origin of mode

*W*is not yet clear. From the character of the doughnut-shape intensity envelope, we suggest that mode

*W*might be a whispering-gallery-like mode that was discussed in Ref. [9]. We also note that the threshold-gain discrimination between the individual modes, otherwise equal to zero in infinitely periodic structures, may result from the differences in their in-plane losses [23]. A closer examination of the field intensities [Fig. 4(b)] distributed at the boundary (dashed circle) of the laser cavity reveals that the lowest-threshold mode

*C*possesses the smallest amount of energy that escapes from the boundary (i.e., the in-plane losses). In contrast, mode

*A*exhibits a much higher threshold because of the relatively larger in-plane losses. This intuitive explanation is confirmed by a quantitative study of the in-plane losses based on the formula derived in Ref. [23]. We find that, for the laser structure specified in Fig. 4, the in-plane losses dominate the threshold and that the in-plane losses of mode

*A*are a factor of 2.7 larger than those of mode

*C*.

#### 3.3. Comparison of experimental and theoretical results and discussion

### 3.3.1. Comparison of band structures and FFPs

Figure 5(a) shows the measured band structure of a fabricated device with the same structural parameters as specified above. The band structure was mapped out around the Γ point by measuring the angle-dependent spontaneous emission spectra well below the lasing threshold [26]. The threshold current (*I*_{th}) at room temperature continuous wave (CW) operation was 25 mA. The measurements were performed in the Γ–*X* and Γ–*J* directions under a CW current level of 0.9*I*_{th} and in the measurements the sample was attached on a heat sink. The dashed curves in Fig. 5(a) replot the data shown in Fig. 3(a) in order to enable a direct comparison with the experimental results. It can be seen that the measured and calculated band structures are in satisfactory agreement. The slight discrepancy may be mainly attributed to the experimental details, such as the change in refractive index owing to carrier injection and thermal effects. We note that in the measured band structure, several additional bands (e.g., bands existing between modes *D*_{1} and *C* near the Γ point) were also observed. These bands are considered to be TM modes as discussed in Ref. [24]. As the laser structure described in Table 1 is not completely symmetric in the vertical direction, TM mode usually coexists with TE mode and may derive some optical gain from the multiple quantum-well active layer. However, since the optical gain is greater for the TE mode than for the TM mode [4], the existence of the TM bands rarely becomes an issue for our laser devices. It is also interesting to note that, especially near the Γ point, the intensity of spontaneous emission from the antisymmetric modes *A*–*C* is extremely weak compared to that from modes *D*_{1} and *D*_{2}. This is due to the difference in their radiative nature, as presented in Fig. 3(b). The lasing spectrum at the Γ point measured at 1.2*I*_{th} current level is shown in Fig. 5(b), which evidently suggests that single-mode lasing indeed occurs at band-edge mode *C* (yellow dashed line). This observation is consistent with the lasing mode predicted above.

Next, we examine the FFP and polarization profile of the lasing mode. (See Ref. [23] for details of the calculation method for FFP.) Figure 6(a) shows the calculated FFP of mode *C*. A doughnut-shape beam with a divergence angle of around 1° is obtained; this closely matches the measured FFP shown in Fig. 6(b). The polarization profiles shown in the insets indicate that mode *C* exhibits an azimuthal polarization, which stems from the monopolar nature of its field pattern, as depicted in Fig. 3(c).

### 3.3.2. Lasing mode control by tuning air-hole size

As has been demonstrated in the previous works on triangular-lattice PC-SELs [3, 12], single-mode lasing action may also be achieved at band-edge mode *A* (i.e., the hexapole mode). Based on our CWT analysis, we find that by tuning some structural parameters, e.g., by varying the air-hole size, mode *A* may be favored for lasing. As an example, we list in Table 3 the normalized threshold gain of the low-threshold modes *A*–*C* for air holes with a larger filling factor of *f* = 0.26. The data given in Table 3 suggest that mode *A* is the lasing mode with a threshold-gain discrimination of around 20 cm^{−1} against the competing modes (i.e., modes *B*_{1} and *B*_{2}). The calculated FFP and polarization profile of mode *A* are shown in Fig. 7(a). It is seen that mode *A* exhibits a doughnut-shape pattern with hexagonal symmetry and possesses a characteristic polarization: The electric-field vector is rotated 4*π* around the circumference of the doughnut beam (instead of 2*π* for the monopole mode *C*). The weak side lobes of the calculated FFP result from the radiation field being terminated abruptly in the numerical calculation.

To further demonstrate the predictive power of our CWT analysis, we fabricated another device with *f* = 0.26. By measuring the band structure and the lasing spectrum in the same manner as described for Fig. 5, we confirmed that lasing indeed occurs at band-edge mode *A*. The measured FFP and polarization profile are shown in Fig. 7(b), and they are in good agreement with the calculated FFP shown in Fig. 7(a). The slight asymmetry in the measured FFP (e.g., the weak streak in the central part of the doughnut-shape beam) can be attributed to the nonuniform electrical pumping or the asymmetric geometry of the gain region determined by carrier diffusion in practical experimental situations.

## 4. Conclusions

A full 3D CWT was developed for triangular-lattice PC-SELs with TE polarization by extending our previous CWT model [24]. When compared to the previous model, the key feature of the extended theory is that six basic waves were introduced in the coupled-wave equation to describe the *C*_{6v} symmetry of the eigenmodes in the vicinity of the second-order Γ point; additionally, both surface radiation and in-plane losses, which are critical for analyzing the finite-size laser structure, were included. By solving the extended coupled-wave equation, modal properties, including the band structure, radiation constant, threshold gain, field intensity profile, and FFP were studied. Comparison of the calculated band structure and FFP of the lasing mode with experimental results showed good agreement. Single-mode lasing at the monopole mode *C* predicted by our threshold gain analysis was confirmed by experimental data. Furthermore, it was shown that the hexapole mode *A* can be favored for lasing by using a larger air-hole size; this was also confirmed by the experiments.

Although the analyses performed in this work were restricted to circular air holes, the power of our theory lies in the fact that it can be applied to air holes of any arbitrary shape as it inherently incorporates a large number of high-order waves [21]. Moreover, the general coupled-wave equation derived in this work is applicable not only to triangular lattice but also to centered-rectangular lattices with *C*_{2v} or *C*_{4v} symmetry, thereby unifying the analysis of the modal properties of PC-SELs. Our future work will focus on a more systematic study of the modal properties of PC-SELs by varying structural parameters of PCs, such as the air-hole shape and the lattice angle, to obtain a comprehensive understanding the physics of PC-SELs and exploit their potential applications.

The presented 3D-CWT analysis is a *linear* theory, and it gives reliable solutions only at or near the threshold current. When the current level is far above the threshold, modal properties of the practical laser devices may deviate from predictions based on the linear analysis. We will leave the above-threshold 3D-CWT analysis that incorporates the optical nonlinearities such as hole burning, nonuniform gain distribution, and thermal effects to future work.

## Appendix A: Solutions of radiative and high-order waves for triangular-lattice PCs

The radiative and high-order waves [i.e., *E _{x,m′,n′}* and

*E*in Eq. (15)] can be solved in a manner similar to that described in our previous works [21,24]. The major difference is that the number of basic waves that serve as sources to excite radiative and high-order waves are six in this case, instead of four. Therefore, we shall only show the final analytical solutions.

_{y,m′,n′}## Radiative waves

The radiative waves (for which *m* = *n* = 0) can be expressed as

*G*(

*z,z′*) = −

*i*/2

*β*·

_{z}*e*

^{−iβz|z−z′|}is an approximated Green’s function [21,22],

*β*=

_{z}*k*

_{0}

*n*

_{0}(

*z*), and

**V**= (

*R*

_{1},

*S*

_{1},

*R*

_{2},

*S*

_{2},

*R*

_{3},

*S*

_{3})

*.*

^{t}## High-order waves

The integrals of high-order waves (for which |*m*^{2} + *n*^{2}| > 1) included in the second term of the right-hand side of the coupled-wave Eq. (15) can be expressed as

As a consequence, the overlap integrals appearing on the right-hand side of Eq. (15) can be replaced by terms only associated with basic waves. Finally, a coupled-wave equation for infinite periodic structures [for which the derivative terms in Eq. (15) are neglected] can be written in the matrix form as

where**C**is a 6 × 6 matrix. The matrix elements of

**C**can be written as where

*m′*,

*n′*| ≤ 10 [21]. Note that the matrix

**dk**

_{0,mn}represents the variation of the in-plane wavenumbers resulting from deviation of the individual basic wave vectors from the second-order Γ point, and

**C**

*,*

_{b}**C**

*, and*

_{r}**C**

*correspond to the coupling effects of basic, radiative, and high-order waves, respectively.*

_{h}## Appendix B: Staggered-grid finite-difference method on hexagonal grids

The coupled-wave equation for finite systems can be solved by using a staggered-grid finite-difference method. We rewrite the coupled-wave Eq. (15) for the band-edge modes (for which Δ* _{x}* = Δ

*= 0) in the following form:*

_{y}**C**is a 6 × 6 matrix composed of coupling coefficients defined by Eqs. (A7)–(A11).

To model a triangular-lattice PC having a sixfold symmetry, we discretize the computational domain on a hexagonal grid (see, e.g., Sec. 3.7.2 of Ref. [29]). The location of the six basic-wave components and the coupling coefficients in the vicinity of the (*j,k*)th hexagonal cell of the grid is illustrated in Fig. 8(a). Then, Eq. (B1) becomes

*h*is the distance between two adjacent grid points and the coupling coefficients

*C*(1 ≤

_{mn}*m,n*≤ 6) are assumed to be constant throughout the laser cavity. Figure 8(b) shows a schematic of the circular-shape computational domain (yellow shaded region) considered in this work. The computational domain physically corresponds to the gain region determined by the electrode geometry. Because the shape of the electrode at the bottom of the fabricated device is circular (see Ref. [13] for details), we consider a circular computational domain. We define the radius of the circular-shape gain region as

*L*=

*Nh*(

*N*: integer). Schematic depicted in Fig. 8(b) corresponds to the case when

*N*= 2. By multiplying 2

*L*on both sides of Eq. (B2), we obtain

*δ*+

*iα*)

*L*. In this work, we assume the following absorbing boundary condition: The field amplitude of the basic waves always starts from zero at the boundaries. This generalizes a similar concept for the 1D distributed feedback laser structure originally proposed by Kogelnik and Shank [16]. The boundary consists of the set of points lying immediately outside but closest to one half of the circular domain on the side opposite to the wave’s propagation. As an example, the boundary of

*R*

_{3}is indicated by black squares in Fig. 8(b). We investigated the impact of discretization on the solutions of Eq. (B3) by varying

*N*and confirmed that well-converged solutions can be obtained at

*N*= 7 with a short calculation time (∼1 minute with a personal computer).

## Acknowledgments

The authors are grateful to Dr. T. Asano, Dr. A. Oskooi, Dr. Y. Kurosaka, T. Nakamura, and T. Okino for their helpful discussions and valuable suggestions. This work was partly supported by the Core Research for Evolution Science and Technology (CREST) and Consortium for Photon Science and Technology (C-PhoST) programs of the Japan Science and Technology Agency. One of the authors (Y. Liang) is supported by research fellowships from the Japan Society for the Promotion of Science (JSPS).

## References and links

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

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

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

**4. **I. Vurgaftman and J. R. Meyer, “Design optimization for high-brightness surface-emitting photonic-crystal distributed-feedback lasers,” IEEE J. Quantum Electron. **39**, 689–700 (2003). [CrossRef]

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

**6. **E. Miyai, K. Sakai, T. Okano, W. Kunishi, D. Ohnishi, and S. Noda, “Lasers producing tailored beams,” Nature (London) **441**, 946 (2006). [CrossRef]

**7. **H. Matsubara, S. Yoshimoto, H. Saito, Y. Jianglin, Y. Tanaka, and S. Noda, “GaN photonic-crystal surface-emitting laser at blue-violet wavelengths,” Science **319**, 445–447 (2008). [CrossRef]

**8. **M. Kim, C. S. Kim, W. W. Bewley, J. R. Lindle, C. L. Canedy, I. Vurgaftman, and J. R. Meyer, “Surface-emitting photonic-crystal distributed-feedback laser for the midinfrared,” Appl. Phys. Lett. **88**, 191105 (2006). [CrossRef]

**9. **Y. Chassagneux, R. Colombelli, W. Maineult, S. Barbieri, H. E. Beere, D. A. Ritchie, S. P. Khanna, E. H. Linfield, and A. G. Davies, “Electrically pumped photonic-crystal terahertz lasers controlled by boundary conditions,” Nature (London) **457**, 174–178 (2009). [CrossRef]

**10. **L. Mahler and A. Tredicucci, “Photonic engineering of surface-emitting terahertz quantum cascade lasers,” Laser Photonics Rev. **5**, 647–658 (2011).

**11. **Y. Kurosaka, S. Iwahashi, Y. Liang, K. Sakai, E. Miyai, W. Kunishi, D. Ohnishi, and S. Noda, “On-chip beam-steering photonic-crystal lasers,” Nat. Photonics **4**, 447–450 (2010). [CrossRef]

**12. **S. Iwahashi, Y. Kurosaka, K. Sakai, K. Kitamura, N. Takayama, and S. Noda, “Higher-order vector beams produced by photonic-crystal lasers,” Opt. Express **19**, 11963–11968 (2011). [CrossRef] [PubMed]

**13. **S. Iwahashi, K. Sakai, Y. Kurosaka, and S. Noda, “Centered-rectangular lattice photonic-crystal surface-emitting lasers,” Phys. Rev. B **85**, 035304 (2012). [CrossRef]

**14. **H. Y. Ryu, M. Notomi, and Y. H. Lee, “Finite-difference time-domain investigation of band-edge resonant modes in finite-size two-dimensional photonic crystal slab,” Phys. Rev. B **68**, 045209 (2003). [CrossRef]

**15. **M. Yokoyama and S. Noda, “Finite-difference time-domain simulation of two-dimensional photonic crystal surface-emitting laser,” Opt. Express **13**, 2869–2880 (2005). [CrossRef] [PubMed]

**16. **H. Kogelnik and C. V. Shank, “Coupled-wave theory of distributed feedback lasers,” J. Appl. Phys. **43**, 2327–2335 (1972). [CrossRef]

**17. **M. Toda, “Proposed cross grating single-mode DFB laser,” IEEE J. Quantum Electron. **28**, 1653–1662, (1992). [CrossRef]

**18. **K. Sakai, E. Miyai, and S. Noda, “Coupled-wave model for square-lattice two-dimensional photonic crystal with transverse-electric-like mode,” Appl. Phys. Lett. **89**, 021101 (2006). [CrossRef]

**19. **K. Sakai, J. Yue, and S. Noda, “Coupled-wave model for triangular-lattice photonic crystal with transverse electric polarization,” Opt. Express **16**, 6033–6040 (2008). [CrossRef] [PubMed]

**20. **K. Sakai, E. Miyai, and S. Noda, “Coupled-wave theory for square-lattice photonic crystal lasers with TE polarization,” IEEE J. Quantum Electron. **46**, 788–795 (2010). [CrossRef]

**21. **Y. Liang, C. Peng, K. Sakai, S. Iwahashi, and S. Noda, “Three-dimensional coupled-wave model for square-lattice photonic-crystal lasers with transverse electric polarization: A general approach,” Phys. Rev. B **84**, 195119 (2011). [CrossRef]

**22. **C. Peng, Y. Liang, K. Sakai, S. Iwahashi, and S. Noda, “Coupled-wave analysis for photonic-crystal surface-emitting lasers on air-holes with arbitrary sidewalls,” Opt. Express **19**, 24672–24686 (2011). [CrossRef] [PubMed]

**23. **Y. Liang, C. Peng, K. Sakai, S. Iwahashi, and S. Noda, “Three-dimensional coupled-wave analysis for square-lattice photonic-crystal lasers with transverse electric polarization: Finite-size effects,” Opt. Express **20**, 15945–15961 (2012). [CrossRef] [PubMed]

**24. **C. Peng, Y. Liang, K. Sakai, S. Iwahashi, and S. Noda, “Three-dimensional coupled-wave theory analysis of centered-rectangular lattice photonic crystal with transverse-electric-like mode,” Phys. Rev. B **86**, 035108 (2012). [CrossRef]

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

**26. **K. Sakai, E. Miyai, T. Sakaguchi, D. Ohnishi, T. Okano, and S. Noda, “Lasing band-edge identification for a surface-emitting photonic crystal laser,” IEEE J. Sel. Areas Commun. **23**, 1335–1340 (2005). [CrossRef]

**27. **M. Koba and P. Szczepanski, “The threshold mode structure analysis of the two-dimensional photonic crystal lasers,” Prog. Electromagn. Res. **125**, 365–389 (2012). [CrossRef]

**28. **K. Forberich, M. Diem, J. Crewett, U. Lemmer, A. Gombert, and K. Busch, “Lasing action in two-dimensional organic photonic crystal lasers with hexagonal symmetry,” Appl. Phys. B **82**, 539–541 (2006). [CrossRef]

**29. **A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, 3rd ed. (Artech House, Norwood, 2005).