## Abstract

We present a theoretical analysis of lasing action in photonic crystal surface-emitting lasers (PCSELs). The semiclassical laser equations for such structures are simulated with three different theoretical techniques: exact finite-difference time-domain calculations, an steady-state ab-initio laser theory and a semi-analytical coupled-mode formalism. Our simulations show that, for an exemplary four-level gain model, the excitation of dark Fano resonances featuring arbitrarily large quality factors can lead to a significant reduction of the lasing threshold of PCSELs with respect to conventional vertical-cavity surface-emitting lasers. Our calculations also suggest that at the onset of lasing action, most of the laser power generated by finite-size PCSELs is emitted in the photonic crystal plane rather than the vertical direction. In addition to their fundamental interest, these findings may affect further engineering of active devices based on photonic crystal slabs.

© 2011 Optical Society of America

## 1. Introduction

Ever since the demonstration of the first laser fifty years ago [1], controlling the spontaneous and stimulated emission of atoms, molecules or electron-hole pairs, by placing them in the proximity of complex dielectric or metallic structures has been seen as an attractive way to tailor the lasing properties of active media [2–5]. With the rapid advance of micro- and nano-fabrication techniques over the past decade, along with the development of novel theoretical techniques, this approach has led to the emergence of a number of new coherent light sources, with dimensions and emission properties much different from those found in conventional lasers [6–32].

Among these systems, photonic crystals (PhCs) slabs [33, 37] have already demonstrated their significant potential to enable efficient integrated laser sources [8, 10, 12, 16–19, 22, 23, 26, 31]. In this context, the unique properties of PhCs to achieve simultaneous spectral and spatial electromagnetic (EM) mode engineering can be exploited in two different ways. On one hand, one can take advantage of the extremely high *Q*/*V*_{mode} microphotonic cavities (*Q* and *V*_{mode} being the corresponding quality factor and the modal volume, respectively) that can be introduced in PhCs simply by inducing local variations in the geometry or the dielectric constant of an otherwise perfectly periodic structure [33, 35]. Alternatively, one can benefit from the multidirectional distributed feedback effect occurring at frequencies close to the band-edges in a defect-free PhC slab; this enables coherent lasing emission and polarization control over large areas [10]. In this work, we focus on the latter class of structures, usually referred to as photonic-crystal surface-emitting lasers (PCSELs) [26].

One of the most intriguing aspects of periodic PhC slabs is the presence of guided resonances, the so-called Fano resonances [33, 36, 37]. The physical origin of these resonances lies in the coupling between the guided modes supported by the slab and external plane waves, this coupling being assisted by the Bragg diffraction occurring in the considered structures due to the periodic modulation of their dielectric constant. Importantly, in the case of an infinitely perfectly periodic PhC slab, and essentially due to symmetry reasons [38], some of these Fano resonances are completely de-coupled from the external world (i.e. their *Q*-factor diverges despite lying above the light line). As we show below, in actual finite-size PhC slabs, these *dark* Fano resonances retain some of the properties of their infinite periodic counterparts and, thus, in the limit of large size PhC slabs, can display arbitrarily large *Q* values. Since these dark modes typically have photon lifetimes much longer than those of other modes, we expect them to dominate the lasing properties of PCSELs.

The purpose of this paper is to analyze theoretically how lasing action in PCSELs is influenced by the presence of the dark Fano resonances, using three different theoretical techniques suitable for analyzing microstructured lasers: a generalized finite-difference time-domain method (FDTD) [9, 19, 40, 41], a steady-state ab-initio laser theory [21, 25, 42], and a semi-analytical coupled-mode theory [33,43]. Each of these techniques offers different kinds of insight into the systems. Our results suggest that the physical origin of the low lasing thresholds observed in actual PCSELs, compared to conventional vertical-cavity-surface-emitting lasers (VCSELs) [46], can be explained in terms of dark Fano states. Likewise, we find that, for the exemplary PCSELs structures discussed in this paper, at input pump rates close to the threshold, the presence of dark Fano modes leads to most of the laser power to be emitted in the plane of periodicity rather than in the vertical direction. In comparison with previous work in this area [10, 12, 22, 26, 31], the findings reported in this paper expand the current understanding of lasing action in PCSELs, which, to our knowledge, has hitherto been explained only in terms of the low-group velocity band-edge modes supported by this class of lasers.

The paper is organized as follows. Section 2 discusses the computational methods used in this paper. In this section, we first introduce the general semiclassical approach we have employed to simulate lasing action in PCSELs. Then, we provide a brief summary of each of the three numerical methods used to solve the semiclassical laser equations. Section 3 presents a detailed analysis of the physics of lasing action in different classes of PCSEL structures. Finally, in Section 4 we provide a set of conclusions of this work.

## 2. Methods

#### 2.1. General framework

Before proceeding with the description of the different electrodynamic techniques used in this work, we first present the general theoretical framework from which these techniques are developed.

We model a dispersive Lorentz active medium using a generic four-level atomic system [9, 19, 34]. The population density in each level is given by *N _{i}* (

*i*= 0, 1, 2,3). Maxwell’s equations for isotropic media are given by

*∂*

**B**(

**r**,

*t*)/

*∂t*= −∇ ×

**E**(

**r**,

*t*) and

*∂*

**D**(

**r**,

*t*)/

*∂t*= ∇ ×

**H**(

**r**,

*t*), where

**B**(

**r**,

*t*) =

*μμ*

_{o}**H**(

**r**,

*t*),

**D**(

**r**,

*t*) =

*ɛɛ*

_{o}**E**(

**r**,

*t*) +

**P**(

**r**,

*t*) and

**P**(

**r**,

*t*) is the dispersive electric polarization density that corresponds to the transitions between two atomic levels,

*N*

_{1}and

*N*

_{2}. The vector

**P**introduces gain in Maxwell’s equation and its time evolution can be shown [47] to follow that of a homogeneously broadened Lorentzian oscillator driven by the coupling between the population inversion and the external electric field. Thus,

**P**obeys the equation of motion

*stands for the linewidth of the atomic transitions at*

_{m}*ω*= (

_{m}*E*

_{2}–

*E*

_{1})/

*h̄*, and accounts for both the nonradiative energy decay rate, as well as dephasing processes that arise from incoherently driven polarizations.

*E*

_{1}and

*E*

_{2}correspond to the energies of

*N*

_{1}and

*N*

_{2}, respectively.

*σ*is the coupling strength of

_{m}**P**to the external electric field and Δ

*N*(

**r**,

*t*) =

*N*

_{2}(

**r**,

*t*) –

*N*

_{1}(

**r**,

*t*) is the population inversion driving

**P**. Positive inversion is attained when Δ

*N*(

**r**,

*t*) > 0, in which case the medium is amplifying; when Δ

*N*(

**r**,

*t*) < 0, the medium is absorbing. In order to model realistic gain media, only conditions favorable to the former are considered.

The atomic population densities obey the rate equations

*N*(

**r**,

*t*) > 0, the plus sign corresponds to radiation while the minus sign represents excitation.

*R*

_{p}is the external pumping rate that transfers electrons from the ground state to the third excited level, and is proportional to the incident pump power.

*τ*is the nonradiative decay lifetimes from level

_{ij}*i*to

*j*(

*i*>

*j*) so that the energy associated with the decay term $\frac{{N}_{i}}{{\tau}_{ij}}$ is considered to be lost. Since the total electron density ${N}_{\text{tot}}={\sum}_{i=0}^{3}{N}_{i}\left(\mathbf{r},t\right)$ is conserved in the rate equations, the active medium modeled this way becomes saturable at high values of

*R*

_{p}.

Lasing action in this class of active media is obtained as follows. Electrons are pumped from the ground-state level *N*_{0} to *N*_{3} at a rate *R*_{p}. These electrons then decay nonradiately into *N*_{2} after a short lifetime *τ*_{32}. By enforcing *τ*_{21} ≫ (*τ*_{32}, *τ*_{10}), a metastable state is formed at *N*_{2} favoring a positive population inversion between *N*_{2} and *N*_{1} (i.e. Δ*N* > 0), which are separated by energy *h̄ω _{m}*. In this regime, a net decay of electrons to

*N*

_{1}occurs through stimulated emission and nonradiative relaxation. Lastly, electrons decay nonradiatively and quickly to

*N*

_{0}. Lasing occurs for pumping rates beyond a given threshold ${R}_{\text{p}}^{\text{th}}$, once sufficient inversion (gain) is attained to overcome the total losses in the structure.

In reality, some of these *N _{i}* levels may stand for clusters of closely spaced but distinct levels, where the relaxation processes among them are much faster than that with all other levels. More generally, although we will focus on the particular case of four-level gain media, we expect the lasing properties of the PCSELs to be influenced primarily by the EM properties of the passive dielectric structure, rather than the microscopic details of the gain mechanism. Hence, our results should be broadly applicable to any active device describable by semiclassical laser theory. The specific effects of optimizing lasing action will depend on the gain medium. Thus, as we show below, for a four-level gain medium, it leads to an arbitrary reduction of the lasing threshold; whereas in a semiconductor laser, it may lead to lasing thresholds close to the limit set by the transparency condition.

#### 2.2. Finite-difference time-domain simulations of active media

The finite-difference time-domain method [40] is used to numerically solve for the optical response of the active material. Similar to the derivation of Yee [48], Maxwell’s equations are approximated by the second order center differencing scheme so that both three dimensional (3D) space and time are discretized, leading to spatial and temporal interleaving of the electromagnetic fields. Pixels of the Yee lattice are chosen to be smaller than the characteristic wavelength of the fields. Moreover, this uniform space grid (i.e. Δ*x* = Δ*y* = Δ*z*, where Δ*x*, Δ*y* and Δ*z* are the space increments in the *x*, *y* and *z* directions) should also be made fine enough so that the PCSEL structures in consideration are well-represented. For numerical stability, Von Neumann analysis places an upper bound on the size of the time step, Δ*t* ≤ *S*Δ*x*/*c*, where *c* is the speed of light. The Courant factor *S* is typically chosen to be
$\frac{1}{2}$ to accommodate 3D simulations.

Following Yee, we evolve the **E** and **H** fields at alternate time steps. For simplicity, we show explicitly the set of discretized equations implemented for a one dimensional (1D) setup assuming a non-magnetic and isotropic medium, and denote any functions of space and time as *F ^{n}*(

*i*) =

*F*(

*i*Δ

*z*,

*n*Δ

*t*).

**H**is first updated as

**P**at

*n*+ 1 from the two previous instances of

**P**, and the previous

*N*and

_{i}**E**according to Eq. (1). Note that components of

**P**reside at the same locations as those of

**E**.

**H**and

**P**values to retrieve

**E**at

*n*+ 1:

*N*at

_{i}*n*+ 1 requires

*N*at

_{i}*n*and both the previous and updated

**E**and

**P**values at

*n*and

*n*+ 1. Since the population densities of the four levels are interdependent in Eq. (2) to (5), they must be solved simultaneously by setting up the following system of equations:

**Ã**and

**B̃**are tensors that couple the population in the four atomic levels and updated

*N*values at

_{i}*n*+ 1 in Eq. (11) can be computed by inverting

**Ã**. Note that the atomic population density

*N*is a scalar and depends generally on all three components of

_{i}**E**and

**P**. The above cycle is repeated at each time step until steady state is reached, allowing the full temporal development of the laser mode to be tracked.

In this FDTD scheme with auxiliary differential equations for **P** and *N _{i}*, all physical quantities of the active material are tracked at all points in the computational domain and at all times. The numerical results computed in this

*ab initio*way are exact apart from discretization of space-time, and allow for nonlinear interactions between the media and the fields. In our simulations, the electric, magnetic and polarization fields are initialized to zero except for background noise while the total electron density is initialized to the ground state level. The computational domain is truncated with Bloch periodic boundary conditions or perfectly matched layers (PML), which are artificial absorbing material designed so that the computational grid’s boundaries are reflectionless in the limit Δ

*z*→ 0 [49]. The resolution is consistently set to 20 for every lattice constant,

*a*, where we checked that the relative differences in frequency and

*Q*values between the 20 pixels per

*a*case and the 40 pixels per

*a*case is less than 5%.

#### 2.3. Coupled-mode theory formalism applied to lasing media

Coupled-mode theory (CMT) [33, 43] has been extensively used to study a broad range of different problems in photonics, both in the linear and nonlinear regimes [50–56]. Here, for the first time to our knowledge, we extend that class of analysis to the case of micro-photonic structures that include active media.

Consider first an arbitrary gain medium embedded in an EM cavity with a single lasing mode. The power input to this system, *P*_{in}, can be seen as the result of the rate of work performed by the current induced in the system by the active media, **J**(**r**, *t*), against the electric field of the cavity, **E**(**r**,*t*). Noticing that **J**(**r**, *t*) actually comes from the temporal variation of the polarization density, i.e. **J**(**r**,*t*) = *∂***P**(**r**,*t*)/*∂t* [see definition of **P** in Eq. (1)], *P*_{in} can be written as [57]

Now, we assume that the spatial and temporal dependence in both the electric field and the polarization density can be separated as **E**(**r**, *t*) = **E**_{0}(**r**)*a*(*t*)exp(−*iω _{c}t*) and

**P**(

**r**,

*t*) =

**E**

_{0}(

**r**)

*P*(

*t*)exp(−

*iω*), respectively.

_{m}t**E**

_{0}(

**r**) is the normalized cavity mode profile (∫

*d*

^{3}

*rɛ*

_{0}

*n*

^{2}(

**r**)|

**E**

_{0}(

**r**)|

^{2}= 1) and

*a*(

*t*) is the corresponding slowly-varying wave amplitude, normalized so that |

*a*(

*t*)|

^{2}is the energy stored in the resonant mode [43].

*ω*and

_{c}*ω*stand for the resonant frequency of the cavity and the considered atomic transition [see Eq. (1)], respectively.

_{m}If we further assume that all the power emitted by that system is collected by a waveguide evanescently coupled to the cavity, and that *ω _{c}* =

*ω*, using first order perturbation theory in Maxwell’s equations, one can show from energy conservation arguments [43] that the temporal evolution of the electric field amplitude

_{m}*a*(

*t*) is governed by the equation

*τ*and

_{ex}*τ*are, respectively, the decay rates due to

_{IO}*external*losses (mainly absorption and radiation losses) and due to the decay into the waveguide. The confinement factor ${\xi}_{1}=(1/2){\int}_{A}{d}^{3}r{\left|{\mathbf{E}}_{0}\left(\mathbf{r}\right)\right|}^{2}$, where

*A*denotes the active part of the structure, accounts for the fact that only the active region drives the temporal evolution of the cavity mode amplitude [as seen in Eq. (13)].

We now turn to the analysis of the polarization amplitude *P*(*t*). Within the first-order perturbation theory approach we are describing here, from Eq. (1) one can obtain a simple first-order differential equation for the temporal evolution of *P*(*t*) by making three main assumptions: *(i)* we assume that the linewidth of the atomic transition is much smaller than its frequency, i.e. Γ* _{m}* ≪

*ω*;

_{m}*(ii)*we apply the rotating-wave approximation (RWA), i.e. we consider just the terms that oscillate as exp(±

*iω*); and,

_{m}t*(iii)*we apply the slowly-varying envelope approximation (SVEA) [47], i.e. we assume

*d*

^{2}

*P*(

*t*)/

*d*

^{2}

*t*≪

*ω*(

_{m}dP*t*)/

*dt*. For the structures that we will consider below, these approximations give good agreement between CMT and FDTD.

Thus, the equation of motion for *P*(*t*) is given by

*N*(

*t*)〉 is defined as the population inversion Δ

*N*(

**r**,

*t*) averaged over the mode profile in the gain region of the system

*σ*and

_{m}*ω*are the same as those give in Eq. (1).

_{m}Finally, by projecting Eqs. (2)–(5) onto the cavity mode profile intensity |**E**_{0}(**r**)**|**^{2} over the active region *A*, following a derivation similar to that given above and after some straightforward algebra, one can obtain the following set of equations governing the time evolution of the average population densities 〈*N _{i}*(

*t*)〉 (using the same definition of average as in Eq. (15), and with

*i*=0,...,3)

*ξ*

_{2}is given by ${\xi}_{2}=\left(2{\int}_{A}{d}^{3}r{\left|{\mathbf{E}}_{0}\left(\mathbf{r}\right)\right|}^{4}\right)/\left({\int}_{A}{d}^{3}r{\left|{\mathbf{E}}_{0}\left(\mathbf{r}\right)\right|}^{2}\right)$. Thus, starting from the knowledge of the electric field profile of the resonant cavity

**E**

_{0}(

**r**), and their corresponding decay rates, one can compute all the relevant physical quantities characterizing lasing action in such structure just by solving the linear system of first-order differential equations given by Eqs. (13), (14), (15)–(19). In particular, once such system of equations have been solved, the total emitted power can be easily computed from

*P*(

_{e}*t*) = 2|

*a*(

*t*)|

^{2}/

*τ*. Although the generalization of the approach described here to the case where more than a single mode is lasing is straightforward, for simplicity, we only consider the single-mode case in this paper.

_{IO}The formalism presented here also allows us to explicitly retrieve analytic expressions for the lasing threshold and the slope in the region near threshold, with the spatial contents of the setup entirely embedded in *ξ*_{1} and *ξ*_{2}. For the four-level system considered, assuming *τ*_{10}, *τ*_{32} ≪ *τ*_{21} and that 1/*τ*_{21} ≫ *R*_{p} at steady state, the threshold population inversion is

*τ*

_{tot}is defined so that 1/

*τ*

_{tot}= 1/

*τ*+ 1/

_{IO}*τ*and that all time dependences have been dropped since we are considering the steady state. Correspondingly, the pumping threshold is found to be

_{ex}*(i)*spatially uniform field

*(ii)*only the loss channel related to the cavity’s

*Q*is present, i.e.

*τ*

_{tot}=

*τ*

_{IO}*(iii)*the gain medium fills the whole space considered

*(iv)*

*σ*=

_{m}*ɛ*

_{m}ɛ_{o}λ^{3}

*ω*/4

_{m}*π*

^{2}

*τ*

_{spont}as derived in [47], where

*τ*

_{spont}is the radiative spontaneous lifetime of the lasing transition (i.e. between

*N*

_{1}and

*N*

_{2}in our four-level system).

#### 2.4. Steady-state ab-initio laser theory

Aside from FDTD and CMT, we also employ a newer method, Steady-state Ab-initio Laser Theory (SALT) [21, 25, 42]. In this frequency-domain approach, the laser equations are converted into a set of coupled non-linear wave equations, which are then solved self-consistently to obtain the various multi-mode, steady-state laser properties, including lasing frequencies, thresholds, power output, and the fields inside and outside the cavity. In contrast to CMT, the openness of the cavity is treated exactly instead of using phenomenological loss factors, and the slowly-varying envelope approximation is not employed. SALT was originally formulated to find the stationary solutions of the Maxwell-Bloch equations, which assume a gain medium of two-level atoms, but it can be straightforwardly generalized to treat four-level lasing by appropriate redefinitions of pump and relaxation parameters [44]. In the current work we adapt the four-level formulation of SALT to the classical four-level polarization model of Eqs. (1)–(5); this is the first application of the theory to a novel four-level laser system.

The SALT solves the non-linear lasing equations exactly (except for the use of the rotating wave approximation) in steady-state for single-mode lasing up to and including the second lasing threshold; for multimode lasing, it uses the stationary population approximation (*dN _{i}*/

*dt*= 0) [21], and even in this regime it has been found to agree very accurately with time-domain simulations [45]. A major advantage of the theory is that it achieves higher precision than FDTD with less computational effort as it requires no time-stepping. We will limit our use of the SALT to two dimensional PCSELs, for which the electric and polarization fields are scalars, and will assume single-mode lasing (although we also examine the possibility of a second lasing threshold for completeness). The single-mode SALT begins by re-writing the fields using the periodic ansatz

*ω*is to be determined self-consistently [21, 25, 42], and can include “line-pulling” effects if the gain center

_{L}*ω*differs from the resonance frequency of the passive cavity. In the present work, the gain center will be tuned exactly to the cavity frequency, so we will find that

_{m}*ω*=

_{L}*ω*, but the SALT can also handle the more general detuned case.

_{m}We further assume that the populations are stationary, so the rate Eqs. (2)–(4) reduce to

*ɛ*, and an active part that describes the effect of the inverted gain medium. The nonlinearity arises from the spatial hole-burning term

*h*(

*r⃗*); in the single-mode regime this describes the saturation effect of the mode pattern on the gain medium, reducing the inversion and thus changing its own spatial profile. In the multimode regime, this term would also describe mode coupling, which increases the thresholds of the higher-order modes.

Since Eq. (28) describes laser emission, it must be solved with the outgoing boundary condition, which states that Ψ(**r**) reduces at infinity to a superposition of purely outgoing waves with frequency *ω _{L}*. This non-Hermitian boundary condition is required for a rigorously correct description of an open, steady-state system containing an amplifying medium [21, 42, 44].

The nonlinear problem is solved with the aid of a basis set of “constant flux” (CF) states. For any given *ω*, the CF states {*u _{n}*(

*r⃗*;

*ω*)|

*n*= 1, 2,···} are the discrete solutions to

*η*’s are complex eigenvalues that can be roughly interpreted as various values of the complex dielectric constant which can produce a resonance pole (purely outgoing solution) at the given real

_{n}*ω*[42]. Equation (30) can be solved by variants of existing numerical techniques; in the present paper, we use the finite-element method (FEM). It can be shown that the CF states are self-orthogonal:

*ω*=

*ω*. Hence, the lasing threshold is found by sweeping over a range of frequencies near the gain center

_{L}*ω*, and finding CF eigenvalues obeying

_{m}*δ*

_{0}yields the threshold mode, and the matching

*ω*is the lasing frequency

*ω*at threshold. From ${\delta}_{0}={\delta}_{0}\left({R}_{\text{p}}^{\text{th}}\right)$, we can obtain ${R}_{\text{p}}^{\text{th}}$ using Eq. (27).

_{L}Above threshold, the full SALT expands the lasing mode in the CF basis, accounting for non-linear variations in the mode profile and the lasing frequency. Here, we will instead use a simplification known as the “single-pole approximation (SPA)”, which is valid for high-Q cavities [42]. This theory (SPA-SALT) identifies the lasing mode and frequency with the threshold lasing state, neglecting the *R*_{p}-dependence of *ω _{L}* and the mode profile:

*I*(

*R*

_{p}) is the mode intensity and

*u*is the CF state corresponding to the mode at threshold. SPA-SALT still includes the non-linear saturation effect on the lasing mode, but by eliminating the need to recompute

_{L}*ω*and the lasing mode profile for each

_{L}*R*

_{p}, the above-threshold SALT calculations is greatly speeded up. It gives excellent results for high-Q laser cavities such as PCSELs, although for low-Q cavities (e.g. random lasers) it is much less accurate [42]. Inserting Eq. (33) into Eq. (28), and using Eq. (31), now yields the simple approximate expression

*L*is the out-of-plane height, and we have assumed that the gain medium is distributed uniformly over the region

_{z}*A*. Aside from small differences in the integrands, this expression agrees with the CMT result derived in Eq. (22). In particular, it states that the mode intensities are approximately linear in

*R*

_{p}for ${R}_{\text{p}}>{R}_{\text{p}}^{\text{th}}$.

## 3. Results and discussion

#### 3.1. Passive properties—bandstructures

We begin by studying the dispersion relations of the three PhC lasing structures shown in the three insets of Fig. 1: a 1D cavity structure with 1D periodicity, a 2D slab structure with 1D periodicity, and a 3D slab structure with 2D periodicity (see insets of Figs. 1(a), 1(b), and 1(c), respectively). In all three cases, the corresponding dispersion relations were computed through FDTD calculations by setting up a unit cell of the PhC and imposing PML absorbing boundary conditions on the top and bottom surfaces of the computational domain. Bloch periodic boundary conditions on the electric fields were imposed on the remaining surfaces perpendicular to the slabs. Figure 1(a) illustrates the band diagram of a structure resembling a conventional VCSEL [46], which extends uniformly to infinity in the *x* and *z* directions. In this system, a one-wavelength thick cavity with *n* = 3.55 (e.g. as in InGaAsP) is enclosed by 25 and 30 bilayers of quarter-wave distributed Bragg reflectors (DBRs) on the top and bottom sides of the structure, respectively. The dielectric mirrors consist of alternate layers of dielectric with *n* = 3.17 and *n* = 3.51 (e.g. as in the InP-based lattice matched InP and InGaAlAs, which offers a relatively larger refractive index contrast of Δ*n* = 0.34 at 1.55 *μm* wavelength; allowing broadband, high reflectivity and low penetration depth DBRs to be attained with fewer layers). Pink shaded regions in Fig. 1(a) represent the continuum of bands guided in the DBRs, while the red line represents the air light line that separates the modes that are propagating in air from those that are evanescent in air. Only transverse magnetic (TM) modes with electric field oriented along the *z* direction are considered.

The lowest guided mode [shown as a green line in Fig. 1(a)] is bounded by the light line of the *n* = 3.55 center layer (not shown) and that of the multilayer cladding (bottom edge of the lower continuum region). Thus, this mode is guided within the cavity layer via total internal refraction, just as in regular dielectric waveguide slab, with no means of coupling to air. It is the portion of the second mode which lies above the air light line [plotted as a blue line in Fig. 1(a)] that is useful for laser operation. In fact, it is most often desirable to operate at the frequency that corresponds to *k _{x}* = 0 (the so-called Γ point) so that the power is vertically emitted through the surface in the longitudinal (

*y*) direction. This mode resides in the lowest photonic bandgap of the periodic claddings and, therefore, is trapped within the cavity layer by the high reflectivities (> 96%) of the DBRs. From our calculations, we find that

*Q*, which measures the loss of the VCSEL in the

*y*direction, is 7500 at

*k*= 0 and may generally be increased further by adding more bilayers of the claddings. Thus, VCSEL structures similar to the one described here, resemble a conventional laser cavity in which the eigenmodes are formed in the longitudinal direction due to feedback from the dielectric mirrors and in which the number of the modes increases with the cavity thickness. Notice that the group velocity (

_{x}*v*=

_{g}*dω*/

*dk*) is near zero for small values of

_{x}*k*, which maximizes the wave-matter interaction inside the cavity and, at the same time, enhances the lateral modal confinement.

_{x}Figure 1(b) and 1(c) render the dispersion relations of *air-bridge* type PhC slabs with 1D corrugation and punctured 2D square lattice of air cylinders respectively. These PhC slabs can support Fano resonances. As mentioned in the introduction, these guided resonances appear in the system when periodic air perturbations, introduced in an otherwise uniform slab, enable the coupling between the guided modes supported by the slab and the external radiation continuum, with the strength of this coupling measured by *Q* of the slab structures. One major difference between these PhC slabs and VCSEL-like structures is that in the former light confinement occurs in the in-plane periodic directions due to Bragg diffractions, and in the out-of-plane direction due to index guiding. It is this presence of index guiding in the third dimension that limits the photon lifetime at frequencies above the air light line, leading to far-field radiation. Since discrete translational symmetries exists due to in-plane periodicity, the projected band diagrams are plotted with respected to the lateral wave vectors along the irreducible Brillouin zone. We shall briefly examine the geometries of the two slab structures separately, before drawing the similarities between them when operated as band-edge mode lasers.

The 2D PhC slab sketched in the inset of Fig. 1(b) consists of a 0.3*a*-thick *n* = 3.17 (e.g. as in InP) slab with a set of 1D periodic grooves along the *x*-direction. These grooves are 0.15*a* deep and 0.1*a* wide, and extend uniformly in the *z* direction. Only modes with electric field oriented along *z* are considered. On the other hand, the PhC slab shown in inset of Fig. 1(c) consist of a 0.3*a*-thick *n* = 3.17 (e.g. as in InP) slab punctured with a 2D square lattice of circular air cylinders in the lateral directions, with both depth and radius being equal to 0.25*a*. In this case, only transverse-electric-like (TE-like) modes, with the electric fields primarily horizontal near the center of the slab, are excited. As in the case of the VCSEL, the modes above the light line at Γ are the most desirable for lasing, since they allow the power to be coupled vertically out of the slab surface. Moreover, in this structure, the zero in-plane group velocity facilitates formation of standing waves, as in any conventional cavity, leading to lateral feedback of the eigenmodes. In fact, in the finite size devices that we will be considering next, Δ*k*_{||} ≠ 0 so that the dispersion curves near Γ may be well approximated by the second order Taylor expansion, in which case, *v _{g}* becomes directly proportional to the curvature of the bands. Hence, flat dispersion curves having high density of photonic states and low

*v*are favorable for enhancing light-matter interaction, which is essential for lasing to take place. Note that a VCSEL, on the other hand, has the same direction of periodicity, feedback, and power emission.

_{g}As pointed out in [12], the first set of modes at Γ are ideal for orthogonal out-of-plane surface emission lasers. These are the two band-edge modes shown in Fig. 1(b) and the four band-edge modes shown in Fig. 1(c). The former corresponds to the phase matching condition
${k}_{x}^{\text{d}}={k}_{x}^{\text{i}}+qK$, where
${k}_{x}^{\text{d}}$ and
${k}_{x}^{\text{i}}$ are the diffracted and incident wave vectors respectively, *K* = 2*π*/*a* is the Bragg grating vector, and we only consider *q* = 1 to ensure vertical outcoupling. All other higher lying frequency modes result in additional out-of-plane emission directions at oblique angles from the slab surfaces. For the corrugated slab, the phase matching conditions in the reciprocal space also implies that the waves traveling in the +*x* direction are coupled to those in the −*x* direction within each unit cell, forming an in-plane feedback mechanism, similar to a 1D cavity. These lateral standing waves are in turn coupled into *y* because the Bragg condition is also satisfied along the slab normal, enabling perpendicular surface emission. For the slab shown in Fig. 1(c), phase matching at Γ again couples waves in the four equivalent Γ – *X* directions of a unit cell to the waves emitting in *z*. Here, the main feedback mechanism is provided separately by waves traveling in the ±*x* and ±*y* directions. Further coupling of waves between these orthogonal directions is facilitated by higher order waves traveling in the Γ – *M* directions (see inset of Fig. 1(c) for the definitions of directions in the reciprocal space of a square lattice). Due to the ease of fabrication resulting from the connected nature of the defect-free lattice, as well as other advantages mentioned at the beginning of this section, PhC slab structures hold great potential as laser devices. The key is its ability to excite a single lateral and longitudinal mode over a large 2D lasing area, as a result of multidimensional distributed feedback mechanism described above. Intuitively, we may treat each unit cell as an individual cavity in-sync with its neighbors, to produce coherent laser oscillations, and desired properties of the lasing mode may be affected simply by tuning the design of each lattice cell. This approach has been experimentally realized to control the polarization of the lasing mode [10].

In this work, we focus on another property of the PhC slab that allows it to operate as a high *Q*, low threshold laser: the existence of band-edge modes with infinite photon lifetime, i.e. with no means of coupling out of the slab. This phenomenon occurs for the lower band edge in Fig. 1(b) and for singly degenerate modes in 2D periodic PhC slabs, corresponding to the two lowest band-edge modes at Γ in Fig. 1(c). The absence of radiative components at these points in the band diagram is a result of in-phase superpositions of the forward and backward traveling waves, with in-plane electric field vectors adding destructively. This same feature can be explained using the symmetry mismatch existing between the guided modes in the PhC slab and the diffracted radiation field in air [38]. We shall reinforce these arguments in the next section based on the electric field profiles of the radiation components. Infinite *Q* above the air light line can only be achieved in PhC slabs, this property being absent in VCSELs, or conventional microcavity structures that use high reflectivity mirrors for mode trapping.

In order to study the mode trapping capabilities of the slab structures for use as photonic crystal surface-emitting lasers (PCSEL), we first examine in detail the corrugated slab design. This 2D design, though analytically and computationally much less demanding, encompasses the same essential physics as a 3D PhC slab realizable in experiments, which we will also study at the end of this section.

Figure 2(a) presents the *Q* of the two bands above the light line at small values of *k _{x}*, plotted against frequency, in the vicinity of the bandgap for the PCSEL structure shown in Fig. 1(b), with grooves 0.1

*a*wide and 0.15

*a*deep. We see from the figure that the two band-edge modes differ drastically. The

*Q*of the lower frequency mode diverges rapidly as

*k*→ 0, while that of the next-ordered band remains finite. This is clearly illustrated by the electric field profiles in the unit cell, depicted in the two leftmost panels of Fig. 2(c) for the lower (left) and upper (right) band-edge modes. The unbounded

_{x}*Q*mode, whose radiative electric field component is anti-symmetric about the groove, interferes destructively with itself in the far-field, resulting in no net outcoupling to air. For

*k*away from Γ, this symmetry mismatch is lost, and

_{x}*Q*decreases rapidly but remains large. On the other hand, the second mode is symmetric and vertical emission out of the slab is possible. Note that despite this leakage, most of the electric field is confined within the slab, forming a standing wave pattern due to the lateral feedback mechanism described previously, a signature of Fano resonances. Apart from mode symmetries, the resonances in the slabs are also influenced by the size of the grooves, which may be regarded as periodic dielectric perturbations to an otherwise uniform slab. Results for 1D periodic grooves with depth 0.1

*a*and 0.05

*a*are also shown in Fig. 2(a). Consistent with predictions from the perturbation theory [33], the bandgap decreases with the grooves size while

*Q*increases, approaching the slab waveguide limit of infinity when no grooves are present.

#### 3.2. Passive properties—finite structures

The symmetry of the lower band-edge mode, which forbids outcoupling, is exact only for the infinite (periodic) structure. In any finite system, the photon lifetime is large but finite. Figure 2(b) shows the *Q* factor, as a function of frequency, for finite slabs with lateral sizes ranging from 20 to 320 periods. These results were obtained from FDTD calculations, with the boundary of the computational domain padded with absorbing boundary conditions (PMLs) to mimic the behavior of a slab in free space. A couple of key observations are in order: *(i)* The lower band-edge mode of the finite PhC slab no longer possesses an unbounded *Q*, owing to the fact that an additional loss channel is opened up: energy can now leak from the sides of the slab. This can be observed in the top right panel of Fig. 2(c) for a 20*a* long PhC slab. These lateral losses dominate in the lower band-edge mode. The bottom right panel of Fig. 2(c) shows the symmetric mode, where both vertical and lateral power emission appears equally dominant. It is thus no surprise that the net *Q* of the lower band-edge mode remains higher than that of the symmetric mode [see Fig. 2(b)]. *(ii)* The *Q* of the lasing structure increases with the number of periods, so the lasing threshold correspondingly decreases. We shall quantify the losses in Fig. 3, as functions of the number of periods. *(iii)* The resonant frequencies of the upper band-edge mode are different in the finite and infinite slabs, due to the presence of lateral losses in the former. In the finite system, increasing *ω* leads to a corresponding increase in *k _{x}* and

*v*, and hence a decrease in

_{g}*Q*. In the infinite system, there are no lateral losses, so

*Q*increases with frequency near the band-edge. For the lower band-edge, mode symmetry considerations ensure that

*Q*remains a maximum for both infinite and finite slabs.

Next, we quantify the *Q* values of the corrugated slab in order to understand how the lateral size of the device, *L _{x}*, affects the outcoupling of Fano resonances. Figure 3 compiles the total

*Q*(

*Q*

_{tot}) of the two band-edge modes presented in Fig. 2(b) for PCSEL structures having 0.05

*a*deep grooves, with

*L*ranging from 20

_{x}*a*to 320

*a*. In order to operate the device at typical optical communication wavelength (∼ 1.55

*μ*m), we set

*a*= 675 nm here and in subsequent results. Since a larger PhC slab provides longer confinement time,

*Q*increases with

*L*for both symmetric and anti-symmetric modes. The anti-symmetric mode has higher

_{x}*Q*, due to its reduced vertical emission, as already observed in Fig. 2(c). For

*L*> 100

_{x}*μ*m, the total

*Q*of both modes tends towards that of their infinite counterpart [see Fig. 2(b)]:

*Q*

_{sym}saturates at 1964, whereas

*Q*

_{anti-sym}is unbounded. Therefore, the anti-symmetric mode holds great potential for low-threshold laser operation. Using approximate analytic relationships that govern

*Q*’s dependence on

*L*(unique for each mode), curves fitted to the calculated data are also plotted in Fig. 3. We shall specify these relationships in the next paragraph.

_{x}To better explore the potential of the PhC slab as a vertical emission laser, we decompose *Q*_{tot} into two *Q* values governing the in-plane (*Q*_{||}) and orthogonal out-of-plane (*Q*_{⊥}) directions. The former is a measure of lateral losses from the sides of the slabs while the latter indicates the degree of vertical emissions. They are related by 1/*Q*_{tot} = 1/*Q*_{||} + 1/*Q*_{⊥}. From *Q*_{||} = *ωτ*_{||}/2 (*τ*_{||}/2 is the photon lifetime before escaping from the sides), it can be shown that near Γ, *Q*_{||} scales approximately with the finite slab’s size as
${C}_{1}{L}_{x}^{2}$, where *C*_{1} is a constant independent of *L _{x}*. This scaling may be derived by first quadratically approximating the band near the band-edge as
$\omega \propto {k}_{||}^{2}$ so that the

*v*=

_{g}*dω*/

*dk*

_{||}∝

*k*

_{||}. In addition, taking the limit at Δ

*k*

_{||}Δ

*x*=

*C*, where

*C*is a constant and Δ

*x*=

*L*here, it may be concluded that Δ

_{x}*k*

_{||}scales as 1/

*L*. This sets the characteristic scale for

_{x}*k*

_{||}and hence,

*v*∝ 1/

_{g}*L*. Thus, together with the distance

_{x}*L*the photon travels,

_{x}*τ*

_{||}scales as ${L}_{x}^{2}$. Since both band-edge modes possess low group velocity and thus relatively large lateral photon confinement, and experience the same structural interfaces,

*Q*

_{||}behaves in the same manner for both modes. The same does not apply to

*Q*

_{⊥}, where modal symmetry mismatch considerations act to impede outcoupling. To obtain an approximate scaling of

*Q*

_{tot}with

*L*, we assume ${Q}_{\perp}^{\text{anti-sym}}$ to remain much larger than ${Q}_{\left|\right|}^{\text{anti-sym}}$ while ${Q}_{\perp}^{\text{sym}}$ to be a finite value independent of

_{x}*L*but related to the groove size. Hence, ${Q}_{\text{tot}}^{\text{anti-sym}}\sim {C}_{2}{L}_{x}^{2}$ and ${Q}_{\text{tot}}^{\text{sym}}\sim {Q}_{\text{periodic}}^{\text{sym}}/\left(1+{C}_{3}/{L}_{x}^{2}\right)$, where ${Q}_{\text{periodic}}^{\text{sym}}$ is the value for the symmetric mode of the corresponding infinite slab.

_{x}*C*

_{2}and

*C*

_{3}are constants independent of

*L*. Curve fitting results shown for

_{x}*Q*

_{tot}are made using these relationships and a reasonably good match is achieved. It now becomes clear that great variance of the photon lifetimes for the two modes in Fig. 3 arises purely from their

*Q*

_{⊥}. We conclude that, for small

*L*, the PhC slab behaves as an in-plane emitting device; to excite enough of the band-edge effects and achieve a vertical out-of-plane emitter with large lasing area, it is critical that the dimensions containing the periodicity be made sufficiently large.

_{x}#### 3.3. Lasing—infinite periodic structures

Thus far, we have based our analysis on the properties of the passive dielectric structure. Let us now consider the effects of adding gain; specifically, the four-level gain medium described in Section 2. For simplicity, we assume that the gain is uniformly distributed within the dielectric (the effects of non-uniform gain are outside the scope of this paper, but the present numerical techniques can treat it effectively). For the decay lifetimes, we take *τ*_{10} = *τ*_{32} = 5 ×10^{−14} s and *τ*_{21} = 5 × 10^{−12} s (so a metastable state can form at *N*_{2}). For the coupling constant, we take *σ _{m}* = 1 ×10

^{−4}C

^{2}/kg (this value was obtained assuming that the Purcell effect is negligible); for the total electron density, we take

*N*

_{tot}= 5×10

^{23}m

^{−3}. These values are realistic, and similar to those used in Ref. [9]. For each structure, we choose a different value of the gain center

*ω*, in order to select the mode that we wish to lase; in the FDTD calculations, this frequency is set to the frequency of the corresponding passive mode. The gain linewidth Γ

_{m}*is taken to be 0.002 (2*

_{m}*πc*/

*a*), which is sufficiently narrow to avoid exciting neighboring modes.

First, we compute the lasing properties of the *infinite* slab. The computational domain is similar to the one shown in Fig. 1(b), with periodic boundary conditions along the left and right boundaries and PML absorbers along the top and bottom boundaries. Figure 4 shows the resulting plot of output power versus *R*_{p}. Three different structures, with groove depths of 0.05*a*, 0.1*a*, and 0.15*a*, are simulated; the groove width is kept at 0.1*a*, and slab thickness at 0.3*a*. The filled circles in this plot are the results of FDTD calculations (time-stepping until steady-state laser operation was achieved); the solid lines are the CMT predictions, with parameters fitted from separate FDTD calculations of the passive structure’s Fano resonance frequency, electric field mode profile, and *Q*.

For each of the calculations in Fig. 4, the gain center *ω _{m}* is situated at the resonance frequency of the symmetric mode, as determined by the passive-structure FDTD calculations presented earlier. (Since this is the infinite structure, the anti-symmetric mode does not allow power emission.) For groove depths [0.05

*a*,0.1

*a*, 0.15

*a*], we have

*ω*= [0.435, 0.444, 0.456] ×(2

_{m}*πc*/

*a*) and

*Q*values [1964, 451, 230]; see Fig. 2(a). As expected, the laser threshold is inversely proportional to

*Q*; physically speaking, higher input pump rates are needed to overcome larger losses. Moreover, the three structures exhibit very similar rates of growth of output power,

*dP*/

*dR*

_{p}. As we shall see, this is not true for finite structures.

The agreement between FDTD and CMT is very good. In particular, CMT predicts that the output power grows linearly with *R*_{p} above the lasing threshold, and the FDTD results are very close to linear. The match remains excellent for *R*_{p} as much as an order of magnitude above the lasing threshold. This shows that the CMT model that we have developed greatly complements the FDTD approach. The CMT is particularly useful for *R*_{p} near threshold, where FDTD computations are very time-consuming due to the temporal turn-on delay before lasing action begins. For larger *R*_{p}, the results begin to deviate; the influence of the gain media on the fields can no longer be taken to be linear, so second order corrections to CMT are required and the lasing modes are no longer accurately described by the modes of the linear (passive) cavity.

#### 3.4. Lasing—finite structures

Figure 5 shows *R*_{p} versus total power output obtained by FDTD and CMT, for finite slabs. (In the FDTD calculations, the PML absorber is now placed along all four boundaries of the computational domain.) Three different slab widths are used: *L _{x}* = 20

*a*, 40

*a*, and 80

*a*. We fix the groove depth at 0.05

*a*, with all other parameters kept the same as in Fig. 4.

The left-hand plot in Fig. 5 shows the upper band-edge modes, while the right-hand plot shows the anti-symmetric lower band-edge. Readers are referred to Fig. 2(b) or the figure caption for the *Q* values and frequencies of these modes. As expected, the structures with larger *L _{x}* have lower lasing thresholds, due to stronger diffraction of the waves, which is needed for better feedback and modal confinement. The band-edge effects in finite PhC slabs depend upon the degree of overlap between the lasing modes and the periodic dielectric, as well as the lateral sizes in their periodic planes [26]. The close proximity of the gain medium to the air grooves in our setup ensures that the former condition is well met, so that band-edge mode laser operation can be achieved for

*L*of as little as 20

_{x}*a*. It is also noted from the two plots that the lower band-edge modes give rise to lower thresholds than the upper band-edge modes, for slabs of equal

*L*. Moreover, the slopes

_{x}*dP*/

*dR*

_{p}are different. The large 80

*a*slab, with the largest lasing area, emits the most power and therefore exhibits the highest slope. This can clearly be seen from the insets of Fig. 5(a) and 5(b), where linear plots of the same data are shown for values near the thresholds.

The FDTD calculations were performed using a resolution of 20 pixels per *a*, which is relatively low; in particular, the groove depth of 0.05*a* corresponds to one pixel. The low resolution is due to computational limitations, especially near threshold, where very many time steps are required to bring the laser to its steady state. The CMT parameters were fitted from passive FDTD calculations performed with the same spatial resolution, and the fact that the CMT results agree well with FDTD indicates that the two methods face consistent discretization errors, as expected.

In order to confirm the FDTD/CMT findings, we have also carried out more precise calculations using the SALT described in Section 2.4. Since the SALT is a frequency-domain method, it does not face the time-stepping problem of FDTD near threshold, and it is possible to apply much finer spatial discretization. The CF basis functions (and hence the lasing modes) were computed via an FEM technique, using a non-uniform triangular mesh with maximum mesh size 0.01*a* within the dielectric structure. Convergence was tested by halving the mesh size, whereupon ∼1% deviation in the computed CF eigenvalues was observed.

Figure 6 shows the CF eigenvalue spectrum {*η _{n}*(

*ω*)} of the finite PCSEL structure, for several different values of the total slab length

*L*. Recall from Section 2.4 that the CF eigenvalues are the discrete complex contributions to the dielectric function from the gain medium, required to produce a resonance at a given real frequency

_{x}*ω*. For each

*ω*, we find that one particular CF state has

*η*(

_{n}*ω*) lying significantly closer to the real axis (i.e. requiring less amplification) than all the others. The spatial structure of this CF state, shown in the inset of Fig. 6 for

*L*= 20

_{x}*a*, corresponds closely to the anti-symmetric Fano mode of the passive structure (see Fig. 2). The complex CF eigenvalues depend on

*ω*, and increasing

*ω*causes them to drift leftwards in the complex

*η*plane (reflecting increasing mode confinement). The CF eigenvalue spectra in Fig. 6 are plotted at the optimal threshold lasing frequency

*ω*of each structure. In the SALT calculations, we did not choose

_{L}*ω*using the FDTD calculations of the passive structure, as we did for the FDTD and CMT calculations; this would be inappropriate, as the SALT calculation uses a different and finer grid (this discretization mismatch is further exacerbated by the fact that Γ

_{m}*is chosen to be small). Instead, the thresholds were found using the self-consistent procedure described in Section 2.4, and the optimal choice of gain center*

_{m}*ω*is the one that minimizes ${R}_{\text{p}}^{\text{th}}$ for the anti-symmetric lasing mode. As we see from (32) and in Fig. 6, this corresponds to making the relevant CF eigenvalue

_{m}*η*(

_{n}*ω*) purely imaginary, so that

_{L}*ω*≃

_{m}*ω*.

_{L}Figure 7 shows the effect of the structure’s lateral size *L _{x}* on the optimal lasing frequency

*ω*, the lasing threshold ${R}_{\text{p}}^{\text{th}}$, and the power slope

_{L}*dP*/

*dR*

_{p}. The increase of

*ω*with

_{L}*L*agrees with the FDTD results shown in Fig. 2, while the numerical values of

_{x}*ω*differ from the FDTD predictions by ∼ 4%, a very acceptable deviation considering the difference in resolution between the two calculations. Likewise, ${R}_{\text{p}}^{\text{th}}$ differs from the FDTD result by ∼ 9%, but shows a similar decrease with

_{L}*L*. The power slope, given by Eq. (36), was found to increase approximately linearly with

_{x}*L*. The power output calculated by the SALT is shown by the dashed curves in Fig. 5. Apart from the aforementioned 9% difference in ${R}_{\text{p}}^{\text{th}}$, these results are in good agreement with FDTD, and in particular are in slightly better agreement than CMT for large

_{x}*R*

_{p}.

A unique advantage of SALT is its ability to predict higher modal thresholds. Within SPA-SALT, predicting the second modal threshold, including mode competition, is a simple extension of the single-mode theory [42]. The results of such a calculation, shown as the inset to Fig. 7, indicate that a second lasing mode will never turn on for this choice of *ω _{m}* and Γ

*. For experimental systems, in which the gain parameters are not fully controllable, such calculations are useful for estimating the range of single-mode operation.*

_{m}#### 3.5. 3D slab lasers

We now turn our attention to the 3D slab system illustrated in the inset Fig. 1(c), consisting of a square lattice of air cylinders. The slab thickness is 0.3*a*, while the radius and depth of the air cylinders are 0.15*a*. This design operates based on the same principles as the simpler 1D periodic grooves design, so that the physical concepts explored previously may be equally applied in this case. In Fig. 8, the magnetic and electric field profiles of TE-like excitations are provided for the four modes at Γ, two being non-degenerate (two lowest frequency modes) and the remaining pair is degenerate. Analogous to the anti-symmetric modes that exist in corrugated slabs, the non-degenerate modes of such infinite periodic slabs have infinite photon lifetime, which again may be attributed to mode symmetry mismatch with the radiative continuum [37]: at Γ of the square lattice, its irreducible representation may either be 1D (singly-degenerate) or 2D (doubly-degenerate), and possesses the full symmetry of the lattice. Such symmetries enforced upon the singly-degenerate modes ensure that their in-plane radiative components cancel. This can be seen from the electric field mode profiles in Fig. 8(a) and 8(b), where the *E _{x}* components are of opposite directions relative to the air cylinder. The same holds for

*E*. Hence, no coupling to air is observed. On the other hand, vertical radiation of the electric field occurs for the degenerate modes in Fig. 8(c) indicating low finite

_{y}*Q*values. We apply the CMT approach to calculate the power output from a unit cell for the degenerate mode at three radii of the air cylinders shown in Fig. 8(d). As in the corrugated slab,

*Q*increases for smaller air cylinders leading to lower threshold pump rate. Practical considerations with regards to size of the structural periodic perturbations include the ease of fabrication as dimensions scale down, and also the need for close proximity to the gain layers for enhancement of the band-edge effects.

Lastly, we examine the finite size 3D slab in Fig. 9. Having already verified the CMT predictions with FDTD method, the former is again adopted to calculate the power output of finite 0.3*a*-thick PhC slabs with air cylinders 0.15*a* deep and diameter 0.3*a*. Three sizes of the PhC region are studied: 15*a* × 15*a*, 25*a* × 25*a* and 35*a* × 35*a*. In order to excite PhC states and to model realistic conditions in similar PhC lasers operating with an optically or electrically pumped central area, the 2D PhC region has to be surrounded by un-pumped regions. This is achieved in our simulations by truncating the finite size PhC in air and extending the uniform dielectric slab into the PML, which surround the whole computational domain. The lasing mode considered is that of the first singly-degenerate mode shown in Fig. 8(a) and may be compared to the field profiles presented in Fig. 9 for the slab with 15*a* × 15*a* PhC region. The slope [as defined in Eq. (22)] and threshold pump rate improves for the larger PhC, consistent with what we would have expected, while the PCSEL remains single-mode. The primary losses for the small-sized PhC considered here is through the lateral leakage into the absorbing boundaries. Further note that magnitudes of the output power, and hence slope, are significantly higher than those found for the corrugated slabs in Fig. 5. This can be understood as a consequence of an additional dimension present in the current calculations. We confirm that the singly degenerate modes for these finite slabs have high *Q*s relative to the degenerate modes and hence, should be the ones within the first set of frequencies at Γ that, in practice, are selected for lasing. Such mode selection (singly-degenerate) should be possible owing to the distinct frequencies and field profiles of the four modes at Γ, and its lower pumping rate requirements. We also note here that upon successful coupling to the desired laser modes, *Q* may further be increased by enhancing the confinement geometrically: *(i)* in the lateral directions by employing PhC heterostructures, and *(ii)* in the out-of-plane direction via the addition of DBRs on top or below the PhC slab [58].

## 4. Conclusions

We have shown that the lasing action in PCSELs originates from Fano resonances. Of particular interest for lasing are dark Fano resonances. We used three different theoretical techniques suitable for studying these systems. We have seen that in actual finite-size structures, the large *Q* factors displayed by these dark states lead to a significant reduction of the corresponding lasing threshold with respect to conventional VCSELs. In addition, our calculations suggest that, for input pump rates close to the threshold, PCSEL structures emit most of their lasing power in the plane of periodicity rather than in the vertical direction. However, notice that this lasing power can be directed into the out-of-plane direction simply by perturbing the symmetry of dark Fano resonances [10]. We believe the findings reported here provide further physical insight into lasing action in PCSELs and will help designing active devices based on this class of systems. Finally, it should also be mentioned that spontaneous emission has not been included in our simulations as it does not affect the physics explored and the conclusions of this work. This leads to no output power detected below the threshold pump rate
${R}_{\text{p}}^{\text{th}}$. Spontaneous emission may be effectively simulated in FDTD by introducing noise-like dipole sources in our gain medium [9].

## Acknowledgments

We acknowledge helpful discussions with Bo Zhen, Dr. Peter Bermel, Dr. Ling Lu, Dr. Ofer Shapira, Dr. Ardavan Oskooi, Prof. Steven Johnson, Prof. John Joannopoulos, Prof. Erich Ippen, Prof. Leslie Kolodziejski, and Prof. Hui Cao. This work was supported in part by the MRSEC Program of the National Science Foundation under award number DMR-0819762, the MIT S3TEC Energy Research Frontier Center of the Department of Energy under Grant No. DE-SC0001299, and the Army Research Office through the Institute for Soldier Nanotechnologies under Contract No. W911NF-07-D0004. J.B.-A. acknowledges financial support by MICINN through Ramón y Cajal program, grant no. RyC-2009-05489. Song-Liang Chua acknowledges the financial support from the DSO National Laboratories, Singapore.

## References and links

**1. **T. H. Maiman, “Stimulated optical radiation in ruby,” Nature **187**, 493–494 (1960). [CrossRef]

**2. **E. Yablonovitch, “Inhibited spontaneous emission in solid-state physics and electronics,” Phys. Rev. Lett. **58**, 2059–2062 (1987). [CrossRef] [PubMed]

**3. **S. Haroche and D. Kleppner, “Cavity quantum electrodynamics,” Phys. Today **42**, 24–30 (1989). [CrossRef]

**4. **Y. Yamamoto (ed.), *Coherence, Amplification, and Quantum effects in Semiconductor Lasers* (John Wiley & Sons, 1991), and references therein.

**5. **H. Yokoyama, “Physics and device applications of optical microcavities,” Science **256**, 66–70 (1992). [CrossRef] [PubMed]

**6. **N. M. Lawandy, R. M. Balachandran, A. S. L. Gomes, and E. Sauvain, “Laser action in strongly scattering media,” Nature **368**, 436–438 (1994). [CrossRef]

**7. **D. S. Wiersma and A. Lagendijk, “Light diffusion with gain and random lasers,” Phys. Rev. E **54**, 4256–4265 (1996). [CrossRef]

**8. **O. Painter, R. K. Lee, A. Scherer, A. Yariv, J. D. O’Brien, P. D. Dapkus, and I. Kim, “Two-dimensional photonic band-gap defect mode laser,” Science **284**, 1819–1821 (1999). [CrossRef] [PubMed]

**9. **X. Jiang and C. M. Soukoulis, “Time dependent theory for random lasers,” Phys. Rev. Lett. **85**, 70–73 (2000). [CrossRef] [PubMed]

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

**11. **H. Cao, Y. Ling, J. Y. Xu, C. Q. Cao, and P. Kumar, “Photon statistics of random lasers with resonant feedback,” Phys. Rev. Lett. **86**, 4524–4527 (2001). [CrossRef] [PubMed]

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

**13. **J. C. Johnson, H.- J. Choi, K. P. Knutsen, R. D. Schaller, P. Yang, and R. J. Saykally, “Single gallium nitride nanowire lasers,” Nat. Mater. **1**, 106–110 (2002). [CrossRef]

**14. **D. J. Bergman and M. L. Stockman, “Surface Plasmon Amplification by Stimulated Emission of Radiation: Quantum Generation of Coherent Surface Plasmons in Nanosystems,” Phys. Rev. Lett. **90**, 027402 (2003). [CrossRef] [PubMed]

**15. **X. Duan, Y. Huang, R. Agarwal, and C. M. Lieber, “Single-nanowire electrically driven lasers,” Nature **421**, 241–245 (2003). [CrossRef] [PubMed]

**16. **H.- G. Park, S.- H. Kim, S.- H. Kwon, Y. -G. Ju, J.- K. Yang, J.- H. Baek, S. -B. Kim, and Y.- H. Lee, “Electrically driven single-cell photonic crystal laser,” Science **305**, 1444–1447 (2004). [CrossRef] [PubMed]

**17. **T. Baba, D. Sano, K. Nozaki, K. Inoshita, Y. Kuroki, and F. Koyama, “Observation of fast spontaneous emission decay in GaInAsP photonic crystal point defect nanocavity at room temperature,” Appl. Phys. Lett. **85**, 3889–3891 (2004). [CrossRef]

**18. **H. Y. Ryu, M. Notomi, E. Kuramochi, and T. Segawa, “Large spontaneous emission factor (> 0.1) in the photonic crystal monopole-mode laser,” Appl. Phys. Lett. **84**, 1067–1069 (2004). [CrossRef]

**19. **P. Bermel, E. Lidorikis, Y. Fink, and J. D. Joannopoulos, “Active materials embedded in photonic crystals and coupled to electromagnetic radiation,” Phys. Rev. B **73**, 165125 (2006). [CrossRef]

**20. **C. Gmachl, F. Capasso, E.E. Narimanov, J.U. Nöckel, A. D. Stone, J. Faist, D. Sivco, and A. Cho, “High power directional emission from lasers with chaotic resonators,” Science **280**, 1556–64 (1998). [CrossRef] [PubMed]

**21. **H. E. Türeci, A. D. Stone, and B. Collier, “Self-consistent multimode lasing theory for complex or random lasing media,” Phys. Rev. A **74**, 043822 (2006). [CrossRef]

**22. **B. Bakir, C. Seassal, X. Letartre, P. Viktorovitch, M. Zussy, L. Cioccio, and J. Fedeli, “Surface-emitting micro-laser combining two-dimensional photonic crystal membrane and vertical Bragg mirror,” Appl. Phys. Lett. **88**, 081113 (2006). [CrossRef]

**23. **H. Altug, D. Englund, and J. Vuckovic, “Ultra-fast photonic-crystal nanolasers,” Nat. Phys. **2**, 485–488 (2006). [CrossRef]

**24. **M. T. Hill, Y. -S. Oei, B. Smalbrugge, Y. Zhu, T. de Vries, P. J. van Veldhoven, F. W. M. van Otten, T. J. Eijkemans, J. P. Turkiewicz, H. de Waardt, E. J. Geluk, S. -H. Kwon, Y.- H. Lee, R. Notzel, and M. K. Smit, “Small-divergence semiconductor lasers by plasmonic collimation,” Nat. Photonics **1**, 589–594 (2007). [CrossRef]

**25. **H. E. Türeci, L. Ge, S. Rotter, and A. D. Stone, “Strong interactions in multimode random lasers,” Science **320**, 643–646 (2008). [CrossRef] [PubMed]

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

**27. **S. Gottardo, R. Sapienza, P. D. García, A. Blanco, D. S. Wiersma, and C. López, “Resonance-driven random lasing,” Nat. Photonics **2**, 429–432 (2008). [CrossRef]

**28. **N. I. Zheludev, S. L. Prosvirnin, N. Papasimakis, and V. A. Fedotov, “Lasing spaser,” Nat. Photonics **2**, 351–354 (2008). [CrossRef]

**29. **M. A. Noginov, G. Zhu, A. M. Belgrave, R. Bakker, V. M. Shalaev, E. E. Narimanov, S. Stout, E. Herz, T. Suteewong, and U. Wiesner, “Demonstration of a spaser-based nanolaser,” Nature **460**, 1110–1112 (2009). [CrossRef] [PubMed]

**30. **R. F. Oulton, V. J. Sorger, T. Zentgraf, R.-M. Ma, C. Gladden, L. Dai, G. Bartal, and X. Zhang, “Plasmon lasers at deep subwavelength scale,” Nature **461**, 629 (2009). [CrossRef] [PubMed]

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

**32. **M. P. Nezhad, A. Simic, O. Bondarenko, B. Slutsky, A. Mizrahi, L. Feng, V. Lomakin, and Y. Fainman, “Room-temperature subwavelength metallo-dielectric lasers,” Nat. Photonics **4**, 395–399 (2010). [CrossRef]

**33. **J. D. Joannopoulos, S. G. Johnson, R. D. Meade, and J. N. Winn, *Photonic Crystals: Molding the Flow of Light* (Princeton Univ. Press, 2008).

**34. **A. S. Nagra and R. A. York, “FDTD analysis of wave propagation in nonlinear absorbing and gain media,” IEEE Trans. Antennas Propag. **46**, 334–340 (1998). [CrossRef]

**35. **S. Noda, M. Fujita, and T. Asano, “Spontaneous-emission control by photonic crystals and nanocavities,” Nat. Photonics **1**, 449–458 (2007). [CrossRef]

**36. **A. A. Erchak, D. J. Ripin, S. Fan, P. Rakich, J. D. Joannopoulos, E. P. Ippen, G. S. Petrich, and L. A. Kolodziejski, “Enhanced coupling to vertical radiation using a two-dimensional photonic crystal in a semiconductor light-emitting diode,” Appl. Phys. Lett. **78**, 563–565 (2001). [CrossRef]

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

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

**39. **A. Yariv and P. Yeh, *Photonics: Optical Electronics in Modern Communications* (Oxford University Press, New York, NY, 2007).

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

**41. **A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. D. Joannopoulos, and S. G. Johnson, “MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method,” Comput. Phys. Commun. **81**, 687–702 (2010). [CrossRef]

**42. **Li Ge, D. Y. Chong, and A. D. Stone, “Steady-state Ab Initio Laser Theory: Generalizations and Analytic Results,” arxiv:1008.0628.

**43. **H. A. Haus, *Waves and Fields in Optoelectronics* (Prentice-Hall, 1984).

**44. **L. Ge, “Steady-state Ab Initio Laser Theory and its Applications in Random and Complex Media,” Yale PhD thesis, 2010; also A. Cerjan, Y.-D. Chong, L. Ge and A. D. Stone, unpublished.

**45. **L. Ge, R. Tandy, A. D. Stone, and H. E. Türeci, “Quantitative Verification of Ab Initio Self-Consistent Laser Theory,” Opt. Express **16**, 16895 (2008). [CrossRef] [PubMed]

**46. **H. E. Li and K. Iga, *Vertical-cavity surface-emitting laser devices* (Springer, 2003).

**47. **A. E. Siegman, *Lasers* (Univ. Science Books, 1986).

**48. **K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. **14**, 302–307 (1966). [CrossRef]

**49. **S. D. Glauber, “An anisotropic perfectly matched layer absorbing medium for the truncation of FDTD lattices,” IEEE Trans. Antennas Propag. **44**, 1630–1639 (1996). [CrossRef]

**50. **M. Soljac̆ić, M. Ibanescu, S. G. Johnson, Y. Fink, and J. D. Joannopoulos, “Optimal bistable switching in non-linear photonic crystals,” Phys. Rev. E **66**, 055601 (2002). [CrossRef]

**51. **W. Suh, Z. Wang, and S. Fan, “Temporal coupled-mode theory and the presence of non-orthogonal modes in lossless multimode cavities,” IEEE J. Quantum Electron. **40**, 1511–1518 (2004). [CrossRef]

**52. **J. Bravo-Abad, A. Rodriguez, P. Bermel, S. G. Johnson, J. D. Joannopoulos, and M. Soljac̆ić, “Enhanced non-linear optics in photonic-crystal microcavities,” Opt. Express **15**, 16161–16176 (2007). [CrossRef] [PubMed]

**53. **R. E. Hamam, M. Ibanescu, E. J. Reed, P. Bermel, S. G. Johnson, E. Ippen, J. D. Joannopoulos, and M. Soljac̆ić, “Purcell effect in nonlinear photonic structures: A coupled mode theory analysis,” Opt. Express **16**, 12523–12537 (2008). [CrossRef] [PubMed]

**54. **J. Bravo-Abad, E. P. Ippen, and M. Soljac̆ić, “Ultrafast photodetection in an all-silicon chip enabled by two-photon absorption,” Appl. Phys. Lett. **94**, 241103 (2009). [CrossRef]

**55. **H. Hashemi, A. W. Rodriguez, J. D. Joannopoulos, M. Soljac̆ić, and S. G. Johnson, “Nonlinear harmonic generation and devices in doubly-resonant Kerr cavities,” Phys. Rev. A **79**, 013812 (2009). [CrossRef]

**56. **J. Bravo-Abad, A. W. Rodriguez, J. D. Joannopoulos, P. T. Rakich, S. G. Johnson, and M. Soljac̆ić, “Efficient low-power terahertz generation via on-chip triply-resonant nonlinear frequency mixing,” Appl. Phys. Lett. **96**, 101110 (2010). [CrossRef]

**57. **J. A. Kong, *Electromagnetic wave theory* (EMW Publishing, 2005).

**58. **T. Lu, C. Kao, H. Kuo, G. Huang, and S. Wang, “CW lasing of current injection blue GaN-based vertical cavity surface emitting laser,” Appl. Phys. Lett. **92**, 141102 (2008). [CrossRef]