Analytical coupled-wave model for photonic crystal surface-emitting quantum cascade lasers

An analytical coupled-wave model is developed for surface-emitting photonic-crystal quantum cascade lasers (PhC-QCLs). This model provides an accurate and efﬁcient analysis of full three-dimensional device structure with large-area cavity size. Various laser properties of interest including the band structure, mode frequency, cavity loss, mode intensity proﬁle, and far ﬁeld pattern (FFP), as well as their dependence on PhC structures and cavity size, are investigated. Comparison with numerical simulations conﬁrms the accuracy and validity of our model. The calculated FFP and polarization proﬁle well explain the previously reported experimental results. In particular, we reveal the possibility of switching the lasing modes and generating single-lobed FFP by properly tuning PhC structures.


Introduction
Quantum cascades lasers (QCLs) [1] are unique semiconductor emitters that cover the entire mid-infrared (MIR) region of the electromagnetic spectrum. During the last two decades, significant advances have been made in many aspects of the device, such as layer processing and laser mounting [2,3], thermal and optical loss management [2], as well as the active region design and optimization [3,4]. Nowadays, QCLs are becoming the most promising light sources for many laser-based applications, e.g., trace gas spectroscopy [5], process control [6] and biological sensing [7]. On the other hand, a number of important MIR applications prefer single mode, high output power and good beam quality, including infrared countermeasures [8], remote sensing, and photo-acoustic spectroscopy [9].
A promising approach to keep single-mode lasing within a very large area is to use a twodimensional (2D) periodic structure, i.e., photonic crystal (PhC) [10,11]. This is because PhC enables flexible control of the cavity loss of lasing-mode candidates and therefore allows better control of the mode selection compared to the conventional DFB lasers. Combining advantages of both PhC and intersubband transition of QCLs, MIR PhC-QCLs have been demonstrated first by Colombelli et al. [12] and have increasingly attracted much attention in recent years [13,14].
Despite these experimental advances, systematic theoretical investigation has yet been lacking, particularly for the full three-dimensional (3D) device structure with large periodicity in the PhC plane. This is because it is very challenging to use the commonly-used brute-force simulation tools such as finite-difference time-domain (FDTD) [15,16] or finite-element methods to simulate this type of large-area lasers.
However, previous works have been limited to the infinitely-periodic PhCs and the accuracy is not good enough, as will be pointed out later. In this work, we present an analytical 3D-CWT to accurately model finite-size surface-emitting PhC-QCLs with TM polarization. With this new model, we are now able to investigate various laser properties of interest, including band structure, mode pattern, intensity profile, mode frequency, cavity loss, and far-field pattern (FFP), as well as the cavity loss dependence on device size. The presented 3D-CWT provides efficient treatment of large-area finite-size PhC-QCLs with negligibly small computational time and resources. The validity and accuracy of our model will be verified by comparison with FDTD simulations and previous experimental results.

3D CWT Model
A schematic of the surface-emitting PhC-QCL investigated in this study is displayed in Fig.   1(a). The active region of QCL is based on InGaAs/AlInAs material system and emits MIR light around the wavelength of 8.5 µm [14]. The square-lattice PhC pillars (red) can be are formed by deep-UV lithography and deep dry-etching of the active region (around 2.5 µm thick). Then, the negative space of the PhC pillars is filled with semi-insulating InP (dark gray) using hydride vapour phase epitaxy (HVPE), similar to the fabrication technique of buriedheterostructure PhC-QCLs [14]. We are interested in TM mode because the optical field within the intersubband QCLs is TM-polarized. A typical band structure for TM mode is shown in Fig. 1(b). At the 2nd-order Γ-point, there exist four band-edge modes, referred to as A, B and E (E represents the doubly-degenerate modes E 1 and E 2 ). that dominate the energy inside the PhC layer; these wavevectors are referred to as basic waves (R x , S x , R y , and S y ). They are coupled to each other via two major coupling mechanism: the direct 2D coupling (κ 2 : blue dashed arrow) and direct 1D coupling (κ 3 : red dashed arrow).
The magnetic field H(r) inside the PhC follows Maxwell's equation: where ε(r) is a periodic function, representing the permittivity distribution of the PhC structure, and k is the wavenumber. Inside the PhC layer, the magnetic field of the TM mode can be written as H(r) = (H x (r), H y (r), 0), and can be expanded following Bloch's the- , m x = m + ∆ x , n y = n + ∆ y , m, n are arbitrary integers and ∆ x , ∆ y are deviation from the Γ point [19]. Similarly, the Fourier transform of 1/ε(r) can be written as: permittivity of the active-region material (background InP), and f represents the filling factor (FF) (i.e., the fraction of the area of a unit cell occupied by the active region). For a simple PhC geometry such as circular shape, typically a truncation order of |m, n| ≤ 8 is sufficient for obtaining well-converged solutions [18]. For simplicity, the sidewall of PhC is assumed to be vertical (solution for tilted case is discussed in [17]).
At the 2nd-order Γ-point [18], the Bloch waves inside the PhC layer can be classified into three groups: basic waves ( √ m 2 + n 2 = 1), radiative wave (m = n = 0) and high-order as R x (x, y), S x (x, y), R y (x, y) and S y (x, y). Figure 1(c) schematically illustrates the four basic waves, and two major coupling mechanism among them: direct 2D coupling κ 2 and direct 1D coupling κ 3 , which are explicitly defined in Appendix B.
By taking into account of the envelope function of the Bloch waves within a finite-size structure [22], the coupled wave equation is extended as follows: where δ and α represent the frequency deviation from guided mode and the cavity loss, respectively. The 2nd term on the left-hand side represents the varying envelope of the basic waves, induced by the finite-size effect. The matrix C comprehensively describes the coupling between the basic waves. The details of derivation are presented in Appendix A. Here, the computational domain is assumed as a square-shaped region within where L is the side length. We consider a non-reflection boundary condition . By solving the coupling equation (2) as an eigenvalue problem, we can obtain the mode frequency and cavity loss of the laser mode, as well as the field pattern, intensity profile, FFP, etc.

Results and Discussions
First, to validate our 3D-CWT, we calculate the infinitely-periodic PhC structure. In our calculation model, the permittivity of the active-region and InP are ε a = 3.342 2 and ε b = 3.0637 2 , respectively. Top cladding and substrate materials are also InP (ε b ). The lattice constant a = 2.7 µm, and the thickness of PhC layer t g = 2.5 µm. We assume that both the top and bottom cladding layers (n-doped InP) are thick enough so that the coupling between the waveguide mode and the lossy surface plasmon mode of the metal-semiconductor interface is negligibly small [31]. The calculated band structure for PhC with f = 0.3 is shown in Fig. 2(a). In addition, mode frequency of the band-edge modes dependence on FFs is shown in Fig. 2(b).
Comparison between 3D-CWT and 3D-FDTD indicates an excellent agreement between these two techniques. It should be mentioned that, in Fig. 2(a) our current 3D-CWT shows a significant improvement in accuracy, compared to the previous TM model [29,30]. results [29,30], a discrepancy of δω/ω 0 ≈ 0.5% already occurs at k = 0.02, where δω is the frequency deviation and ω 0 is the center frequency. In contrast, our current results exhibit a much better agreement (δω/ω 0 < 0.05%) even within −0.2 ≤ k ≤ 0.2. This is because a new iteration technique is introduced [24]. Despite of this, the calculation time is still negligibly small (less than 1 minute). Secondly, we calculate the modal property within a finite-size PhC-QCL. By solving Eq.
(2) with the finite-difference method [22], we can directly obtain the normalized cavity loss αL as a function of the normalized mode frequency δL, i.e., mode spectrum. Figure 3 Fig. 3(b), all of which have only one antinode throughout the whole cavity area [23]. Here, the mode intensity profile is defined as: I(x, y) = |R x | 2 + |R y | 2 + |S x | 2 + |S y | 2 [22]. From Fig. 3(b), we see that mode B is the most well-confined mode.
It has the lowest cavity loss (α B = 6.5 cm −1 ) which is more than 5 cm −1 smaller compared to modes E (α E = 11.5 cm −1 ) and A (α A = 21.9 cm −1 ). With such cavity-loss difference (i.e., mode selection), it is possible to achieve stable single-mode lasing operation at mode B.
Here, it is important to note that, the cavity loss calculated by our 3D-CWT consists of both the vertical radiation loss and the in-plane loss. In contrast, the earlier CWT framework developed for TM mode [32,33,34] is essentially a 2D model and therefore cannot correctly model the inplane couplings and vertical radiation loss (which is the output of surface-emitting PhC-QCL). In addition, the smaller mode frequency gap between the individual modes obtained by 3D-CWT indicates that the effective index contrast within a 3D structure is much weaker compared to the 2D case. We also find that, with the increasing of FFs, the lowest cavity-loss mode (usually being selected as the lasing mode) switches from mode A (FF¡20%) to E (25%¡FF¡30%), then to B (30%¡FF¡55%), and finally back to A (FF¿60%). Interestingly, the 2nd and 3rd mode-switching points roughly correspond to the FFs where κ 2 and κ 3 become zero (see Appendix B).
Next, we examine the far-field pattern (FFP) and polarization profiles of the surface emission. By implementing Fourier transform of the radiation field at the surface of the PhC layer [22], we can obtain the FFP and its polarization feature of the individual band-edge modes. In interesting because normally the most favorable lasing mode (i.e., high-Q ⊥ mode A or B) of circular-shaped PhC gives rise to a doughnut beam. The superior FFP of mode E is important for many laser applications requiring single-lobed and high-beam-quality laser beams.
Finally, we investigate the dependence of mode frequency and cavity loss on device size using our developed 3D-CWT. In our calculations, we change the cavity side-length L from 100 a to 2000 a. Figure 6 shows mode frequency and cavity loss as a function of L for FF= 0.25 to the destructive interference and mode E has a finite Q ⊥ -factor due to its radiative nature [35,18]. We can see that, with the increasing of L, both mode frequency and cavity loss become asymptotic to their counterparts of an infinitely periodic structure. Interestingly, we find that the lowest-cavity-loss mode switches from mode E to A for FF= 0.25, whereas for FF= 0.50 mode B maintains to be the lowest-cavity-loss mode.

Conclusion
In this work, an analytical finite-size 3D-CWT model is developed for surface-emitting PhC-QCLs with TM polarization. This model not only provides an efficient treatment of the largearea PhC laser device, but also offers a physical picture to understand the coupling mechanism inside the PhC-QCL structure.
The 3D-CWT model explicitly describes the electromagnetic fields inside the PhC-QCL by considering the couplings between basic waves, radiative waves, and high-order waves. The finite-size effect is included by taking the variation of field envelope into consideration. We demonstrate that, the extended 3D-CWT provides an accurate and reliable solution to a variety of modal properties, including band structure, mode pattern, intensity profile, mode frequency, cavity loss, and radiative beam patterns, etc. The validity of the 3D-CWT is confirmed by comparison with FDTD simulations and previous experimental results. Furthermore, it predicts a flexible control of cavity loss by PhC structural parameters. In particular, it reveals that the lowest-cavity-loss mode (usually the lasing mode) can be selected by properly tuning the FFs and device length. We also find that mode E is suitable for generating single-lobed far-field beam patterns. We believe this work will greatly facilitate the optimization and design of next generation PhC-QCLs.
In the derivation, the high-order derivatives are neglected due to the slow variation of basic wave envelope. Since high-order Fourier terms ξ 1,1 , ξ 1,−1 , ξ −1,1 and ξ −1,−1 are negligibly small compared to ξ 00 , the corresponding terms are also neglected.
By substituting the guided mode equation: into Eq. (5) for all the four basic waves [22], the finite-size coupled-wave Eq. The coupling coefficients κ 2 and κ 3 [see Fig. 1(c)] are particularly important for constructing the 2D coherent resonance within PhC-QCLs with TM polarization [32]. In our 3D-CWT model, κ 2 and κ 3 are explicitly defined by taking into account of the field confinement in the vertical direction: where β 0 = 2π/a, k 0 = 2π/λ, and Γ g represents the percentage of optical energy concentrated in the PhC layer: The dependence of κ 2 and κ 3 on filling factors is shown in Fig. 7. At FF= 0.29 and 0.58, κ 3 and κ 2 become zero, respectively.

C. In-plane loss dependence on cavity area
To investigate the effect of cavity size on the in-plane loss, we perform 2D-FDTD simulations by varying the cavity side-length L. In the simulations, we consider a square-shaped PhC cavity with a side length of L. The cavity is surrounded by perfectly matched layers (PMLs) with thickness of 3 a. An example of the simulated mode B for L = 50 a and FF=0.50 is shown in Fig. 8(a). Figures 8(b) and (c) show the dependence of the in-plane loss (α ) on cavity area L 2 for FF=0.25 and 0.50, receptively. We see that at FF=0.25 mode E has the lowest in-plane loss, whereas at FF=0.50 mode B has the lowest in-plane loss. In addition, we find that the decrease of α is inversely proportional to L 2 ; this is due to the quadratic dispersion of the band structure nearby the 2nd-order Γ-point band edges [36].