Low-threshold lasing action in photonic crystal slabs enabled by Fano resonances

: 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 ﬁnite-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 signiﬁcant 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 ﬁnite-size PCSELs is emitted in the photonic crystal plane rather than the vertical direction. In addition to their fundamental interest, these ﬁndings may affect further engineering of active devices based on photonic crystal slabs.


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][3][4][5]. With the rapid advance of micro-and nanofabrication 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 .
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 bandedges 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 timedomain 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.

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 ∂ 2 P(r,t) where Γ m stands for the linewidth of the atomic transitions at ω 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. σ m is the coupling strength of 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 where ± 1 hω m E· ∂ P ∂t are the stimulated emission rates. For ∆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. τ i j is the nonradiative decay lifetimes from level i to j (i > j) so that the energy associated with the decay term N i τ i j is considered to be lost. Since the total electron density N tot = Σ 3 i=0 N i (r,t) 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 energȳ 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 th p , 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.

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 elec-tromagnetic 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 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 Next, we update the polarization density P at n+1 from the two previous instances of P, and the previous N i and E according to Eq. (1). Note that components of P reside at the same locations as those of E.
P n+1 We can then use these updated H and P values to retrieve E at n + 1: Lastly, N i at n + 1 requires N i at 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: where A and B are tensors that couple the population in the four atomic levels and updated N i values at n + 1 in Eq. (11) can be computed by inverting A. Note that the atomic population density N i is a scalar and depends generally on all three components of 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 spacetime, 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%.

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][51][52][53][54][55][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ω m t), respectively. 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]. ω c and ω m stand for the resonant frequency of the cavity and the considered atomic transition [see Eq. (1)], respectively.
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 = ω m , 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 a(t) is governed by the equation where τ ex and τ IO are, respectively, the decay rates due to external losses (mainly absorption and radiation losses) and due to the decay into the waveguide. The confinement factor ξ 1 = (1/2) A d 3 r|E 0 (r)| 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ω m t); and, (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 Here ∆N(t) is defined as the population inversion ∆N(r,t) averaged over the mode profile in the gain region of the system Notice also that in Eq. (14), the definitions of σ m and ω m are the same as those give in Eq. (1) . 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) where the parameter ξ 2 is given by . 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 /τ IO . 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. 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 where τ tot is defined so that 1/τ tot = 1/τ IO + 1/τ ex and that all time dependences have been dropped since we are considering the steady state. Correspondingly, the pumping threshold is found to be (21) where N tot = Σ 3 i=0 N i . Equation Eq. (21) may be related to the threshold power by P th = R th ph ω A d 3 r N tot . Finally, we define the slope in this paper as the differential ratio of the emitted power over the pumping rate. Our model predicts that Steady state predictions provided by Eqs. Eq. (20), (21), (22) match with results commonly derived in textbooks [39] for a similar four-level system, in the case where the following are assumed in our present model: (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 ω m /4π 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).

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 timedomain 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 The discrete lasing frequency ω L is to be determined self-consistently [21,25,42], and can include "line-pulling" effects if the gain center ω m 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 ω L = ω m , but the SALT can also handle the more general detuned case.
We further assume that the populations are stationary, so the rate Eqs. (2)-(4) reduce to where α(R p ) = 2 τ 21 τ 10 The approximate equalities in the above equations hold for R p ≪ τ −1 21 (≪ τ −1 10 ). Inserting the ansatz Eq. (23) into Maxwell's equations, the polarization equation (1), and the reduced rate equation (24), and making use of the RWA, we obtain Thus, the lasing mode obeys a nonlinear wave equation (28), with a dielectric function consisting of a passive part, ε, 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 with outgoing boundary conditions. The η n '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 ω [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: Comparing Eq. (30) to Eq. (28), we see that the threshold lasing mode corresponds to a CF state with ω = ω L . Hence, the lasing threshold is found by sweeping over a range of frequencies near the gain center ω m , and finding CF eigenvalues obeying Of these, the solution with smallest δ 0 yields the threshold mode, and the matching ω is the lasing frequency ω L at threshold. From δ 0 = δ 0 (R th p ), we can obtain R th p using Eq. (27). 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: where I(R p ) is the mode intensity and u L 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 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 We can obtain the total power output by inserting the above equations back into the physical electric field Eq. (23), and computing the flux of the Poynting vector at infinity. A brief calculation gives, in the physical limit R p ≪ τ −1 21 , the result Here, L z is the out-of-plane height, and we have assumed that the gain medium is distributed uniformly over the region 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 p > R th p .

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 onewavelength 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 x = 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 (v g = dω/dk x ) is near zero for small values of k x , which maximizes the wave-matter interaction inside the cavity and, at the same time, enhances the lateral modal confinement. 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.3a-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.15a deep and 0.1a 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.3a-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.25a. 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 g 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.
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 d x = k i x + qK, where k d x and k i x 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 defectfree 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.1a wide and 0.15a 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 x → 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 (right) band-edge modes. The unbounded 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 x away from Γ, this symmetry mismatch is lost, and 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.1a and 0.05a 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, ap- Fig. 3. Total Q of the two band-edge modes for the finite PhC slab punctured with 0.05a deep grooves, and having total lateral size, L x , ranging from 20 to 320 unit cells. Green lines are the fitted curves using the relationships described in the text and the horizontal line indicates Q value of the corresponding infinite slab (for the symmetric mode case).
proaching the slab waveguide limit of infinity when no grooves are present.

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 20a 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 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 bandedge 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 g , and hence a decrease in 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.05a deep grooves, with L x ranging from 20a to 320a. 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 x for both symmetric and anti-symmetric modes. The anti-symmetric mode has higher Q, due to its reduced vertical emission, as already observed in Fig. 2(c). For L x > 100 µ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 x (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.
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 2 x , where C 1 is a constant independent of L x . This scaling may be derived by first quadratically approximating the band near the bandedge as ω ∝ 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 x here, it may be concluded that ∆k scales as 1/L x . This sets the characteristic scale for k and hence, v g ∝ 1/L x . Thus, together with the distance L x the photon travels, τ scales as L 2 x . 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 x , we assume Q anti-sym ⊥ to remain much larger than Q anti-sym while Q sym ⊥ to be a finite value independent of L x but related to the groove size. Hence, Q anti-sym tot periodic is the value for the symmetric mode of the corresponding infinite slab. C 2 and C 3 are constants independent of L x . Curve fitting results shown for 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 x , 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.

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 ω m , 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π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 Fig. 4. Output power versus R p relationships of the 2D infinite slab described in Fig. 1(b), for three depths of the air grooves at 0.05a, 0.1a, and 0.15a (width remains at 0.1a), with corresponding Q values 1964, 451, and 230 respectively. Both semi-analytic predictions from CMT (solid lines) and FDTD (filled circles) steady-state calculations are plotted for the upper band-edge mode at k x = 0(2π/a). There is good agreement between the semianalytic and calculated values. The threshold is higher for the lower-Q PhC slab which clearly suggests that higher pumping rates are needed to overcome systems with higher losses.
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.05a, 0.1a, and 0.15a, are simulated; the groove width is kept at 0.1a, and slab thickness at 0.3a. 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 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 slabs. Good agreements between the three methods are observed. Slope of the lines changes with L x (see text) while the right plot confirms that the anti-symmetric mode has the largest Q in the finite system that is available for lasing. 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. 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 = 20a, 40a, and 80a. We fix the groove depth at 0.05a, with all other parameters kept the same as in Fig. 4.

Lasing-finite structures
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 x of as little as 20a. 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 x . Moreover, the slopes dP/dR p are different. The large 80a 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.05a 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.01a 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 x . 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 ω. 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 x = 20a, 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 ω L of each structure. In the SALT calculations, we did not choose ω m 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 th p for the anti-symmetric lasing mode. As we see from (32) and in Fig. 6 Fig. 7. (Left) Lasing frequency ω L and optimal gain center ω m , (Center) the lasing threshold R th p , (Right) and the power slope (dP/dR p )/L z , as a function of lateral size L x , computed using the self-consistent ab-initio laser theory (SALT). The inset shows the modal gain of the first and second modes (a mode turns on when its modal gain reaches unity [25]), indicating that the second mode does not turn on at any R p for this choice of gain medium. making the relevant CF eigenvalue η n (ω L ) purely imaginary, so that ω m ≃ ω L . Figure 7 shows the effect of the structure's lateral size L x on the optimal lasing frequency ω L , the lasing threshold R th p , and the power slope dP/dR p . The increase of ω L with L x agrees with the FDTD results shown in Fig. 2, while the numerical values of ω L differ from the FDTD predictions by ∼ 4%, a very acceptable deviation considering the difference in resolution between the two calculations. Likewise, R th p differs from the FDTD result by ∼ 9%, but shows a similar decrease with L x . The power slope, given by Eq. (36), was found to increase approximately linearly with L x . The power output calculated by the SALT is shown by the dashed curves in Fig. 5. Apart from the aforementioned 9% difference in R th p , these results are in good agreement with FDTD, and in particular are in slightly better agreement than CMT for large 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 Γ m . For experimental systems, in which the gain parameters are not fully controllable, such calculations are useful for estimating the range of single-mode operation.

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.3a, while the radius and depth of the air cylinders are 0.15a. 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 y . 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 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.3a-thick PhC slabs with air cylinders 0.15a deep and diameter 0.3a. Three sizes of the PhC region are studied: 15a × 15a, 25a × 25a and 35a × 35a. 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 15a × 15a 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 Fig. 9. Top: Magnetic and electric field profiles corresponding to the singly-degenerate mode in finite 15a × 15a PhC slab structure described in Fig. 8(a). The PhC slab is outlined in green. Top two panels depict the lateral cuts along the xy-plane with magnetic field pointing into the page and electric field pointing to the left, where positive (negative) values are in red (blue). Bottom panels are cuts along the yz-plane with magnetic field pointing downwards and electric field pointing into the page, where positive (negative) values are in red (blue). These results are only for TE-like modes. E y (not shown) has the same profile as E x , except rotated 90 • about z-axis. Bottom: Output power versus R p relationships retrieved from CMT. Frequencies and Q values of the 15a × 15a, 25a × 25a, and 35a × 35a PCSEL structures are 0.431, 0.433, 0.435(2πc/a), and 64, 178, 385, respectively. 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 Qs 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].

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 th p . Spontaneous emission may be effectively simulated in FDTD by introducing noise-like dipole sources in our gain medium [9].