## Abstract

We develop a coupled-wave model that is capable of treating finite-size square-lattice photonic crystal surface emitting lasers with transverse-electric polarization. Various properties of interest including threshold gain, mode frequency, field intensity envelope within the device, far-field pattern, as well as polarization and divergence angle of the output beam for the band-edge modes are calculated. Theoretical predictions of the lowest threshold mode and the output beam profile are in good agreement with our experimental findings. In particular, we show that, contrary to the infinite periodic case, the finite length of the device significantly affects surface emission and mode selection properties of the laser device.

© 2012 Optical Society of America

## 1. Introduction

Two-dimensional (2D) photonic-crystal surface-emitting lasers (PC-SELs) are becoming increasingly important due to their promising functionality and improved performance compared to conventional semiconductor lasers. A number of successful demonstrations underpinning this promise have already been made [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, enabling high-power, single-mode operation [1, 5]. The output beam of such devices is emitted in the direction normal to the 2D PC plane and has a small beam divergence angle of less than 1° due to the large area of coherent oscillation [5, 6]. Furthermore, both the polarization and pattern of the output beam can be controlled by appropriate design of the PC geometry [2,6]. Recent developments of 2D PC lasers have allowed the lasing wavelength to be extended from the near-infrared regime to the mid-infrared, terahertz, and blue-violet regimes [7–12]. In addition, we have recently demonstrated the operation of a PC-SEL with entirely new functionality: on-chip dynamical control of the emitted beam direction, achieved by using a composite PC structure [13].

Despite the experimental advances that have recently been made in the field of 2D PC-SELs, theoretical studies on these types of lasers have thus far been limited. An important but unresolved issue concerns the mechanisms by which the PC structure determines the output characteristics of the device, thereby limiting progress in optimizing the structural design of devices. Computer simulations based on the 2D plane-wave expansion method (PWEM) [3, 14] or the finite-difference time-domain (FDTD) method [15–17], can provide valuable information about the lasing properties of the PC laser cavity. However, there are some inherent limitations to these computational approaches. The 2D PWEM is only applicable to infinite structures, and the FDTD method requires substantial computational resources in order to model finite structures with realistically large areas. A coupled-wave theory (CWT) approach developed by Kogelnik and Shank [18] can circumvent such limitations. However, extending this approach to accurately model PC-SELs is a rather challenging task, since realistic PC-SELs require three-dimensional (3D) analysis. Very recently, we developed a 3D CWT model that affords an exact analytical treatment of the full 3D structure of typical laser devices and achieves a very accurate and efficient analysis of the surface emission properties. This theory incorporates the key issues in modeling surface-emitting-type PC lasers, i.e. the surface emission in the surface normal direction and the in-plane higher-order coupling effects [22], neither of which was appropriately described in the previous works [4, 19–21].

Nevertheless, in our previous study, we restricted our analyses to *infinite*
periodic structures for simplicity [22]. To predict and improve the performance of the practical
PC-SEL device, however, it is essential to consider a *finite-size*
structure. For example, following the arguments of our previous study, two
antisymmetric band-edge modes of PCs with circular air holes should be excited
simultaneously without emitting a laser beam because both of their radiation
constants (i.e. parameters quantifying the surface radiation loss) are shown to be
zero. This is in contrast to the experiments where a laser with a doughnut-shaped
beam pattern operating stably at one of the two antisymmetric modes was demonstrated
[6, 23, 24].
The physical mechanism of this lasing behavior cannot be explained without
considering the *finite-size* effects of the laser device, as we will
show later on. This is our motivation for the present work, where we develop a 3D
CWT model capable of treating the *finite-size* laser device by
extending our previous work [22]. Our objective in this paper is not to exhaustively quantify the
dependence on the large number of parameters that determine device behavior. Rather,
we focus on clarifying the underlying physical mechanism of the effects caused by
finiteness of the device.

The remainder of this paper is organized as follows. Section 2 describes derivations of the coupled-wave equations for finite systems. Section 3 presents analysis results and discusses the finite-size effects. Section 4 concludes with our findings.

## 2. Coupled-wave equations for finite-size structures

In this section, we will present the derivations of the coupled-wave equations for finite systems. We shall not give a complete description of the formulation, but highlight only the major differences from our previous study. For a more detailed derivation and discussion, please refer to our previous paper [22].

The schematic structure of PC-SELs is shown in Fig. 1(a),
in which the electric field can be expressed as **E**(**r**,
*t*) =
**E**(**r**)*e ^{iωt}*. By
eliminating the magnetic field from Maxwell’s equations, we obtain

*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}(\mathbf{r})\simeq {k}_{0}^{2}{n}^{2}(\mathbf{r})+2i{k}_{0}{n}_{0}(z)\tilde{\alpha}(z)$ [22,25], where

*n*

^{2}(

**r**) is a periodic function

*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. Here, we still focus our analysis on the transverse electric (TE) polarization because the lasing mode has been identified as being a TE mode [23]. For TE polarization,

**E**(

**r**) = (

*E*(

_{x}**r**),

*E*(

_{y}**r**), 0) can be expanded according to Bloch’s theorem

*n*

^{2}(

**r**) can be expanded as

*β*

_{0}= 2

*π*/

*a*,

*a*is the lattice constant,

*m*,

*n*are arbitrary integers, and

*ξ*(

_{m,n}*z*) is the high-order Fourier coefficient term. We note that

*ξ*(

_{m,n}*z*) is zero outside the PC region. For simplicity, we assume that air holes within the PC region have perfectly vertical sidewalls (tilted case is discussed in Ref. [26]) such that ${n}_{0}^{2}(z)$ represents the average dielectric constant of PC and that

*ξ*is independent of

_{m,n}*z*within the PC region. By substituting Eqs. (2) and (3) into Eq. (1) and collecting all terms that are multiplied by the factor

*e*

^{−imβ0x−inβ0y}, we obtain

*E*and

_{x}*E*with respect to

_{y}*x*and

*y*. As defined previously [22], the wavevectors can be classified into three groups according to their in-plane wavenumber, $\sqrt{{m}^{2}+{n}^{2}}{\beta}_{0}$: basic waves ( $\sqrt{{m}^{2}+{n}^{2}}=1$), high-order waves ( $\sqrt{{m}^{2}+{n}^{2}}>1$), and radiative waves (

*m*=

*n*= 0). By assuming a separable form of solutions for the fields, we can solve Eqs. (4)–(6) analytically. In the resonant case at the second-order Γ point shown in Fig. 1(b), the four basic waves can be expressed as [22]

*R*(

_{x}*x*,

*y*) and

*S*(

_{x}*x*,

*y*) represent the amplitudes of basic waves propagating in the +

*x*and −

*x*directions, respectively, and likewise,

*R*(

_{y}*x*,

*y*) and

*S*(

_{y}*x*,

*y*) represent the amplitudes of waves propagating in the +

*y*and −

*y*directions, respectively. These four basic waves are assumed to have identical field profiles in the

*z*-direction, denoted by Θ

_{0}(

*z*), which is the same as the field profile of the fundamental guided mode for a multilayer structure without PCs [25]. We express the wave equation for the fundamental guided mode in terms of Θ

_{0}(

*z*) as

*β*is the propagation constant, which satisfies

*β*≃

*β*

_{0}in the vicinity of the second-order Γ point. In this work, we calculate

*β*and Θ

_{0}(

*z*) in Eq. (11) by employing the transfer matrix method [27] and normalize Θ

_{0}(

*z*) as ${\int}_{-\infty}^{\infty}{\left|{\mathrm{\Theta}}_{0}(z)\right|}^{2}dz=1$.

In order to obtain the coupled-wave equations, Eqs. (7)–(10) are substituted into Eqs. (4)–(6) for (*m*, *n*) = {(1, 0), (−1, 0), (0, 1), (0, −1)}. Here, without loss of generality, we focus on the case for which (*m*, *n*) = (1, 0). We then only need to consider Eqs. (7) and (5). Substitution of Eq. (7) into Eq. (5) gives

*R*is assumed to vary slowly compared to exp(−

_{x}*iβ*

_{0}

*x*) so that its second spatial derivative terms in Eq. (5) can be neglected. Then Eq. (11) is substituted into Eq. (12) to yield

*E*

_{y,m′,n′}on the right-hand side of Eq. (13) represents all the waves that may couple to

*R*, including basic, high-order, and radiative waves as described above. Specifically, we express the radiative waves [for which (

_{x}*m*′,

*n*′) = (0, 0)] as

Finally, we can obtain the coupled-wave equation for (*m*, *n*) = (1, 0) by multiplying Eq. (13) by
${\mathrm{\Theta}}_{0}^{*}(z)$ on both sides and integrating over (−∞, ∞) along the *z* direction. Three more coupled-wave equations for (*m*, *n*) = {(−1, 0), (0, 1), (0, −1)} can be derived in analogous fashion. We write the four coupled-wave equations in the following form:

*ω*

_{0}is the Bragg frequency,

*n*is the effective refractive index for the fundamental guided mode, and is the mode loss given by $\alpha =\frac{{k}_{0}}{{\beta}_{0}}{\int}_{-\infty}^{\infty}{n}_{0}(z)\tilde{\alpha}(z){\left|{\mathrm{\Theta}}_{0}(z)\right|}^{2}dz$.

_{eff}*κ*

_{±2,0}and

*κ*

_{0,±2}are the backward coupling coefficients defined as

*ξ*= 0 outside that range.

_{mn}As the fields of the radiative waves (Δ*E _{x}*(

*z*), Δ

*E*(

_{y}*z*)) and the high-order waves (

*E*(

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

*E*(

_{y,m,n}*z*)) are unknown, we cannot yet evaluate the right-hand sides of Eqs. (15)–(18). In order to determine these fields, Eqs. (4)–(6) must be solved for these waves. Details concerning the solutions of Δ

*E*(

_{x}*z*), Δ

*E*(

_{y}*z*),

*E*(

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

*E*(

_{y,m,n}*z*) are shown in Appendix A. Finally, the coupled-wave Eqs. (15)–(18) can be rewritten in matrix form as

**C**is a 4 × 4 matrix (see Appendix A). By considering a finite laser cavity area with 0 ≤

*x*,

*y*≤

*L*(

*L*: laser cavity length in both

*x*and

*y*directions) and defining appropriate boundary conditions, Eq. (21) can be discretized using the staggered-grid finite-difference method [28]. This creates an eigenvalue problem with eigenvectors ${\left({R}_{x}^{j,k},{S}_{x}^{j,k},{R}_{y}^{j,k},{S}_{y}^{j,k},{R}_{x}^{j+1,k},{S}_{x}^{j+1,k},{R}_{y}^{j+1,k},{S}_{y}^{j+1,k},\dots \right)}^{t}$ and normalized eigenvalues (

*δ*+

*iα*)

*L*.

## 3. Analysis results and discussions

In this section, we will show numerical results of the coupled-wave equations derived above. By solving the coupled-wave Eq. (21), we can obtain various properties of interest, including mode frequency, threshold gain, field intensity envelope profile within the device, far-field patterns (including the polarization profile) of the band-edge modes. The effects of finite size on these properties will be discussed.

The structural parameters of the laser structure [see Fig. 1(a)] to be studied are shown in Table 1. The PC layer is embedded inside the multilayer structure, which is assumed to support only a single guided mode. 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). In this work, we focus our analyses on two typical air-hole shape designs [see the inset of Fig. 1(a)]: circular (CC) and equilateral triangular (ET), both of which have been studied experimentally in our previous works [6, 24]. We consider a boundary condition that is defined as

*a*= 295 nm and the air-hole filling factor

*f*= 0.16, the device length

*L*=70

*μ*m (all these values were inferred from experimental data). Here, it is important to note that

*L*=70

*μ*m approximately corresponds to a 240

*a*× 240

*a*lasing area, which is almost impossible to use 3D FDTD method to compute.

#### 3.1. Threshold gain behavior and field intensity envelope

In finite systems, mode frequency and threshold gain correspond to the real part and imaginary part of the normalized eigenvalue of Eq. (21), (*δ* +*iα*)*L*, respectively. As an example, we show in Fig. 2(i) a plot of the normalized threshold gain (*αL*) as a function of deviation from the Bragg condition (*δL*) for a PC-SEL with ET air holes. Though a large number of modes exist due to the numerical calculation, we identify the four band-edge modes (A, B, C and D) using the techniques described in Ref. [21]. First, we evaluate the field intensity envelope of the individual modes throughout the laser cavity, which modulates the fast-varying Bloch waves in the PC lattice and can be determined by [21]

The normalized threshold gains of modes A–D are shown in Table 2 in which the values for the CC air-hole case were obtained in the same manner. We note that modes C and D exhibit a much higher threshold gain compared with modes A and B, which is attributed to the symmetric nature of their electric-field distributions with respect to the air holes as shown in Fig. 2(iii). As the lasing action occurs at the mode with the lowest threshold gain (loss), Table 2 indicates that mode A is favored for lasing for both CC and ET air holes when the device length *L* = 70 *μ*m, with a large threshold-gain discrimination of over 20 cm^{−1} against mode B. This is consistent with the experimental observations [23, 24] where stable single-mode operations at mode A were demonstrated.

It is important to note that the single-mode lasing behavior at mode A cannot be explained based on the infinite solutions [22]. For an infinite periodic system with the vertical structure shown in Table 1, the radiation constants *α _{inf}* of modes A and B are

*α*

_{inf,A}=

*α*

_{inf,B}= 0 cm

^{−1}for CC and

*α*

_{inf,A}= 12.35 cm

^{−1},

*α*

_{inf,B}= 0.82 cm

^{−1}for ET air holes, respectively. From these results, one might expect that both modes A and B are favored for lasing for the CC air holes because both of their radiation constants (losses) are equal to zero; and for the ET air holes, mode B lases due to its much smaller radiation constant compared with that of mode A. Apparently, these predictions are in contrast to our finite analysis results and experimental findings. We suggest that the in-plane loss of finite structures might have a great impact on the mode selection of PC-SELs, which will be discussed in detail in the following sections.

#### 3.2. Far-field pattern and polarization

The far-field pattern and polarization of the output laser beam is one of the most unique features of 2D PC-SELs. As demonstrated in Refs. [2, 6], both far-field pattern and polarization can be tailored by appropriately engineering the PC geometry. Theoretically, it is possible to calculate the far-field patterns of 2D PC-SELs based on our 3D CWT model since the the radiation field emanating from the surface of the finite-size laser structure is formulated explicitly (see Appendix A). The time-dependent amplitude of the far field *F*_{j=x,y} (subscripts denote the *x* and *y* components, respectively), a function of the deflection angle (*θ _{x}*,

*θ*) from the surface normal, can be found by taking the Fourier transform of the radiation field and multiplying by an obliquity factor (cos

_{y}*θ*+ cos

_{x}*θ*− 1) [30],

_{y}*E*(

_{j}*x*,

*y*,

*z*=

*d*/2) is the complex radiation field just above the PC surface [for further detail, see Eq. (A4) in Appendix A],

_{PC}*d*is the depth of the PC layer, and

_{PC}*z*= 0 is defined at the center of the PC layer. Then, the time-averaged far-field intensity (i.e. beam pattern) is given as

*x*and

*y*components of the far field using Eq. (26), we can directly evaluate the effect of the air holes on the polarization of the output beam.

Figures 3(a) and 3(b) show the calculated far-field patterns of mode A (i.e. the lowest threshold mode) for the two air-hole shapes. A doughnut-shaped profile is obtained for CC, whereas a single-lobed profile is obtained for ET air holes. The insets show the *x* and *y* components of the far field, indicating that CC air holes exhibit an azimuthal polarization, and ET ones exhibit an almost linear polarization in the *y* direction. These semi-analytical results closely match our experimental observations shown in Figs. 3(c) and 3(d) [6, 24]. The beam divergence angle of all of the far-field patterns is around 1°, reflecting the large area of coherent oscillation. The small side lobes of the calculated far-field patterns are caused by the abrupt termination of the radiation field at the edges of the laser cavity.

It is interesting to note that the quantity of surface emission for the symmetric CC air-hole shape is zero when an infinite periodic structure is considered [22], whereas it has a finite value for a finite periodic structure. This difference can be understood by considering the effect of finite device length on the electromagnetic field distribution within the device. For an infinite periodic structure, the light from mode A cannot be coupled to the external system because perfect destructive interference occurs everywhere due to the antisymmetric field distribution pattern with respect to the perfectly symmetric air holes. On the other hand, the mode field is spatially restricted for a finite system, leading to a small shift of the electromagnetic field with respect to the air holes. Figure 4(a) shows the electromagnetic field distribution patterns at (i) the center of the device and toward (ii–vi) the edges of the device. Although the mode field is kept antisymmetric with respect to the air holes in the center region, the field in the regions away from the center is no longer antisymmetric with respect to the air holes, thereby resulting in an imperfect destructive interference.

To examine the resultant interference effect in further detail, we plot in Fig. 4(b) the intensities of basic waves propagating in the ±*x* directions (i.e. |*R _{x}*|

^{2}and |

*S*|

_{x}^{2}) and the radiation field intensity along the line

*y*=

*L*/2. The radiation field intensity can be represented by |

*R*+

_{x}*S*|

_{x}^{2}, which reflects the interference effects of the two counter-propagating basic waves

*R*and

_{x}*S*[see Eq. (A4) in Appendix A and Eq. (B11) in Appendix B by noting that Fourier coefficient terms (

_{x}*ξ*) for CC air holes are real numbers]. From Fig. 4(b), we can clearly see that the interference is indeed perfectly destructive at the center but becomes imperfect toward the edges of the device. The situation for ET can be understood in the same manner except that even at the central region the perfect destructive interference is suppressed because of the asymmetric air-hole effects [31].

_{m,n}The divergence angle of the far-field pattern is closely related to the finite size of the device. When the device is infinitely large, the surface emission is strictly in the surface normal direction. However, the realistic device has a finite length *L*, and therefore allows a wavenumber fluctuation of *δk* ≃ 2*π*/*L* from the Γ point, which was demonstrated by 3D FDTD simulations for finite structures [16, 17]. In fact, this wavenumber fluctuation corresponds to the slightly shifted electromagnetic field distributions discussed above [see Fig. 4(a)]. Due to this fluctuation, the emitted beam deviates from the surface normal direction with a deflection angle *δθ*. The relation between *δθ* and the wavenumber fluctuation can be expressed as

*k*

_{0}= 2

*π*/

*λ*(

*λ*: the resonance wavelength in free space). For example, when

*λ*= 960 nm,

*L*=70

*μ*m,

*δθ*≃ 0.8°, which is a good approximation of the results shown in Fig. 3.

#### 3.3. Device length dependence of threshold gain

The analysis for an infinite system corresponds to the case when *L* → ∞. In this case, threshold gain and mode selection properties are determined only by the surface radiation loss (for simplicity, in this work we do not consider absorption loss). However, for a finite structure, influence of the in-plane loss (i.e. the power escaping from the edges of the laser cavity) has to be taken into account. To clarify the effects of the in-plane loss and highlight the difference between finite and infinite analyses, we investigate the threshold dependence on device length.

We plot in Fig. 5 the normalized mode frequency and threshold gain as a function of the device length *L* for CC and ET air holes. Also included as dashed lines are the infinite solutions (*δ* and *α* are first calculated for infinite periodic structures and then normalized by using (*δ* + *β*_{0})/(*n _{eff}*

*β*

_{0}) [22] and

*αL*, respectively). From Figs. 5(a) and 5(b), we note that mode frequency increases with respect to

*L*and converges quickly to the infinite solutions (dashed lines). The threshold gain, in contrast, shows a quite different behavior: as shown in Figs. 5(c) and 5(d), in both cases, the laser threshold gain (solid curves) initially decreases with increasing

*L*and then begins to be asymptotic to the infinite solutions (dashed lines). This behavior can be intuitively understood as follows. For small values of

*L*, the dominant loss mechanism in the laser is the in-plane loss. In particular, at very short lengths, the optical power exits the gain region at the edges of the laser cavity without being diffracted out of the PC surface, leading to a high threshold gain. However, for larger values of

*L*, the in-plane loss steadily decreases and the surface radiation loss becomes dominant. In fact, infinite solutions can correctly predict the lasing mode (i.e. the lowest threshold mode) of the device only when

*L*is large enough, e.g.

*L*> 400

*μ*m (see Fig. 5). For a relatively small

*L*, e.g.

*L*< 100

*μ*m, the influence of the in-plane loss on threshold behavior should not be neglected.

To evaluate the influence of the in-plane loss quantitatively, we derive a formula (Appendix B) to describe the power flow in 2D PC-SELs by extending the energy conservation theorem developed for 1D DFB lasers [32], which can be expressed as

Here,*P*,

_{stim}*P*, and

_{edge}*P*represent the stimulated emission power, the power escaping from the edges of the laser cavity (i.e. the in-plane loss), and the surface radiation power, respectively. Therefore, the percentage of the in-plane loss power with respect to the total stimulated power is given by Table 3 shows the calculated

_{rad}*P*/

_{edge}*P*of modes A and B for CC and ET air holes with different device lengths. We can clearly see that, for both air-hole shapes, when

_{stim}*L*= 70

*μ*m, the in-plane loss is indeed noticeably larger than the surface radiation loss, indicating that in-plane loss plays a very important role in the mode selections of the laser device. While when

*L*=400

*μ*m, the in-plane loss drastically decreases and its influence on mode selections might be neglected. Here, it is important to note that, in the case of ET air holes, when the device length

*L*is increased from 70

*μ*m to 400

*μ*m, the lasing mode likely switches from mode A to mode B which has a smaller radiation constant [see Fig. 5(d)].

## 4. Conclusion

In conclusion, we have developed a CWT model for finite-size square-lattice PC-SELs with TE polarization by extending our previous 3D CWT framework. Using this model, we calculated various properties of 2D PC-SELs, including threshold gain, mode frequency, field intensity envelope within the device, and the profile of the output beam (far-field pattern and polarization). The theoretical predictions of the lowest threshold mode and the output beam profile showed good agreement with our experimental findings, thereby affirming the validity of our analyses.

Our finite-size analysis results demonstrate that finite device length indeed strongly influences surface emission and mode selection properties of the PC-SEL device. Phenomena that cannot be explained by the infinite CWT analysis are clarified by considering the finite-size effects of the device. We showed that the finite device length may lead to a small shift of the electromagnetic field with respect to the air holes. This slightly shifted field suppresses the otherwise perfect destructive interference and thus enhances surface emission from the regions away from the center of the cavity, which accounts for the finite surface radiation loss emitted from the non-radiative (antisymmetric) mode of the CC air holes. By investigating the device length dependence of threshold gain, we find that, at a very large device length *L* (e.g. *L* > 400 *μ*m), the finite-size analysis results become asymptotic to the infinite solutions, while at a relatively small *L* (e.g. *L* < 100 *μ*m), the predicted lasing mode obtained using the finite-size analysis may be quite different from that based on the infinite structures. This is attributed to the influence of the in-plane loss, which can be noticeably large at short device lengths. The in-plane loss can provide a strong mode selection mechanism with a large threshold-gain discrimination of over 20 cm^{−1}, which allows us to explain the stable single-mode operation observed in our previous experiments.

The insights obtained in this work are essential for understanding the performances and physics of a practical PC-SEL device. Moreover, this theory also forms the basis for more complicated analyses of the properties of the device, including slope efficiency and the nonlinear effects [33] such as gain saturation [32], non-uniform gain spatial distribution, and hole burning. We believe that further theoretical work based on this semi-analytical 3D CWT framework will enable efficient optimization of the structures of 2D PC-SELs for a range of applications as well as facilitating a comprehensive understanding of the physics of PC-SELs.

## Appendix A: Solutions of radiative and high-order waves

We assume, following Ref. [25], that: (1) Only basic waves are important in generating radiative and high-order waves; (2) For radiative and high-order waves, their derivatives in Eqs. (4)–(6) are negligible compared to the others; (3) *α̃* is small and thus may be neglected. Under these assumptions, Eqs. (4)–(6) are reduced to

*E*and

_{x,m,n}*E*can be solved in a similar manner as that described in Ref. [22], thus we shall only show the final solutions of radiative and high-order waves.

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

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

*β*=

_{z}*k*

_{0}

*n*

_{0}(

*z*); It should be noted that Δ

*E*(or Δ

_{x}*E*) is a function of

_{y}*x*and

*y*. As a consequence, the second terms of the right-hand sides of the coupled-wave Eqs. (15)–(18) can be replaced by terms only associated with the basic waves.

## High-order waves

The integrals of high-order waves (for which $\sqrt{{m}^{2}+{n}^{2}}>1$) included in the third terms of the right-hand sides of the coupled-wave Eqs. (15)–(18) can be expressed as

## The elements of matrix **C**

Finally, the matrix **C** in Eq. (21) can be written as

**C**

_{1D},

**C**

*and*

_{rad}**C**

_{2D}correspond to the conventional 1D feedback coupling, out-of-plane coupling via radiative waves, and 2D optical coupling via high-order waves, respectively. It should also be noted that

**C**

_{1D}and

**C**

_{2D}are Hermitian matrices, whereas

**C**

*is not an Hermitian one.*

_{rad}## Appendix B: Derivation of energy conservation theorem for 2D PC-SELs

The coupled-wave Eq. (21) can be rewritten as

*C*(1 ≤

_{ij}*i*,

*j*≤ 4) =

*C*

_{1D,ij}+

*C*

_{rad,ij}+

*C*

_{2D,ij}[see Eqs. (A10)–(A13)]. Multiplying Eq. (B1) by ${R}_{x}^{*}$, Eq. (B2) by ${S}_{x}^{*}$, Eq. (B3) by ${R}_{y}^{*}$, and Eq. (B4) by ${S}_{y}^{*}$, and adding the four equations and their conjugates, we obtain

*C*

_{1D,ij}and

*C*

_{2D,ij}terms have vanished due to their Hermitian property, and only the term ${\kappa}_{v,i}=-\mathit{Im}\left\{\frac{{k}_{0}^{4}}{2{\beta}_{0}}{\int}_{PC}G\left(z,{z}^{\prime}\right){\mathrm{\Theta}}_{0}({z}^{\prime}){\mathrm{\Theta}}_{0}^{*}(z)d{z}^{\prime}dz\right\}$ which is closely associated with the out-of-plane coupling remains [see Eq. (A14)]. In order to model the power flow in the 2D PC-SELs, we integrate Eq. (B5) over the whole laser cavity area, i.e., 0 ≤

*x*,

*y*≤

*L*. Then we have

For a laser structure with boundary condition given by Eq. (22), the field amplitudes of basic waves at the four edges can be expressed as

*P*describes the total stimulated power inside the laser structure,

_{stim}*P*represents the power escaping from the edges of the laser cavity (i.e. the in-plane loss), and

_{edge}*P*represents the radiation power emitted from the device surface [note that this quantity is proportional to the intensity of the radiative waves described by Eq. (A4)], respectively. Equation (B8) states the fact that the power generated inside the laser structure is equal to the sum of the power escaping from the edges of the structure and the surface radiation power.

_{rad}## Acknowledgments

The authors are grateful to Dr. K. Ishizaki, Dr. M. Yamaguchi, Dr. A. Oskooi, and Mr. T. Okino for helpful discussions and valuable suggestions. This work was partly supported by the Core Research for Evolution Science and Technology program of the Japan Science and Technology Agency (CREST-JST) and by the Global COE Program. Three of the authors (Y. Liang, C. Peng, and S. Iwahashi) are supported by research fellowships of 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 **441**, 946 (2006). [CrossRef] [PubMed]

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

**8. **G. Xu, Y. Chassagneux, R. Colombelli, G. Beaudoin, and I. Sagnes, “Polarized single-lobed surface emission in mid-infrared, photonic-crystal, quantum-cascade lasers,” Opt. Lett. **35**, 859 (2010). [CrossRef] [PubMed]

**9. **L. Sirigu, R. Terazzi, M. I. Amanti, M. Giovannini, and J. Faist, “Terahertz quantum cascade lasers based on two-dimensional photonic crystal resonators,” Opt. Express **16**, 5206–5217 (2008). [CrossRef] [PubMed]

**10. **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 **457**, 174–178 (2009). [CrossRef] [PubMed]

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

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

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

**14. **M. Plihal and A. A. Maradudin, “Photonic band structure of two-dimensional systems: the triangular lattice,” Phys. Rev. B **44**, 8565–8571 (1991). [CrossRef]

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

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

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

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

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

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

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

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

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

**24. **W. Kunishi, D. Ohnishi, E. Miyai, K. Sakai, and S. Noda, “High-power single-lobed surface-emitting photonic-crystal laser,” Conference on Lasers and Electro-Optics (CLEO), CMKK1, Long Beach, May, 2006.

**25. **W. Streifer, D. R. Scifres, and R. D. Burnham, “Coupled wave analysis of DFB and DBR lasers,” IEEE J. Quantum Electron. **13**, 134–141 (1977). [CrossRef]

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

**27. **M. J. Bergmann and H. C. Casey, “Optical-field calculations for lossy multiple-layer AlGl-1N/InxG1-xN laser diodes,” J. Appl. Phys. **84**, 1196–1203 (1998). [CrossRef]

**28. **D. H. Sheen, K. Tuncay, C. E. Baag, and P. J. Ortoleva, “Parallel implemantation of a velocity-stress staggered-grid finite-difference method for 2-D poroelastic wave propagation,” Comput. Geosci. **32**, 1182–1191 (2006). [CrossRef]

**29. **E. Miyai and S. Noda, “Phase-shift effect on a two-dimensional surface-emitting photonic-crystal laser,” Appl. Phys. Lett. **86**, 111113 (2005). [CrossRef]

**30. **A. Yariv and P. Yeh, *Photonics: Optical Electronics in Modern Communications*, 6th ed. (Oxford University Press, 2007).

**31. **Fundamentally, radiation fields emitted from the center of the laser cavity have similar properties to those emitted from an infinite periodic structure described in Ref. [22]. Unlike CC air holes, Fourier coefficients (*ξ _{m,n}*) of the dielectric function ε(

**r**) for ET air holes are complex numbers. Therefore, the radiation field intensity is proportional to |

*ξ*

_{−1,0}

*R*+

_{x}*ξ*

_{1,0}

*S*|

_{x}^{2}[see Eq. (A4) in Appendix A and Eq. (B11) in Appendix B]. These complex Fourier coefficient terms multiplied to basic waves may change the phase difference of the waves diffracted vertically, resulting in a suppression of the destructive interference.

**32. **H. A. Haus, “Gain saturation in distributed feedback lasers,” Appl. Opt. **14**, 2650–2652 (1975). [CrossRef] [PubMed]

**33. **S. H. Macomber, “Nonlinear analysis of surface-emitting distributed feedback lasers,” IEEE J. Quantum Electron. **26**, 2065–2074 (1990). [CrossRef]