Open Access
Add to BibSonomyAbstract
We show comparison of four different numerical methods for simulating Photonic-Crystal (PC) VCSELs. We present the theoretical basis behind each method and analyze the differences by studying a benchmark VCSEL structure, where the PC structure penetrates all VCSEL layers, the entire top-mirror DBR, a fraction of the top-mirror DBR or just the VCSEL cavity. The different models are evaluated by comparing the predicted resonance wavelengths and threshold gains for different hole diameters and pitches of the PC. The agreement between the models is relatively good, except for one model, which corresponds to the effective index method. The simulation results elucidate the strength and weaknesses of the analyzed methods; and outline the limits of applicability of the different models.
©2010 Optical Society of America
1. Introduction
In the recent years the development of Photonic-Crystal-based Vertical-Cavity Surface-Emitting Lasers (PC-VCSELs) is getting more and more advanced. These devices have strong potential for leaving laboratories into the market due to their unique properties, which make them a perfect choice for many applications. These properties include a stable single-mode operation [1, 2, 3], high-speed modulation [4, 5, 6] and polarization control [7, 8]. Typical PC-VCSELs consist of a classical VCSEL cavity surrounded by Distributed Bragg Reflectors (DBRs) of high reflectivity, with lateral photonic crystal structure which provides lateral light confinement. In most designs additional oxide aperture is placed inside the cavity, however—contrary to classical VCSELs—its role is limited only to limiting the current spreading. The photonic crystal has a form of cylindrical holes located in various parts of the device. In the simplest case—and therefore the most popular one—the holes are etched in the top DBR. However, there are other possibilities like drilling the whole structure or placing the holes solely in the cavity, which can improve some properties of PC-VCSEL but although constitutes a technological challenge.
An efficient application of PC-VCSELs requires a comprehensive modeling of their optical, electrical and thermal properties. Unfortunately the lack of axial symmetry and complex 3D structure of these devices makes this task difficult. Especially the modeling of optical field is a challenging task, as the analyzed structures show large complexity and very high index contrast. Therefore, a strong scattering is taking place at the edges of the holes and a typical linear polarization (LP) approximation is no longer valid. In consequence, application of popular simplified models is impossible, especially in situations where one needs to determine not only the resonant wavelength of the PC-VCSEL cavity, but also the cavity Q-factor or the gain characteristics of the laser.
In this article we compare four optical models: two vectorial, one scalar and one hybrid scalar-vectorial. We use them to predict the emitted light wavelength and the threshold gain of gallium arsenide PC-VCSEL. The schematical diagram of this laser is shown in Fig. 1 and described in details in section 3. The comparison of the models developed and utilized by three unrelated research groups have been performed within a framework of the European COST action MP0702.
Fig. 1. (a) Schematic structure of the analyzed PC-VCSEL. (b) Top view of photonic crystal arrangement.
The paper is organized as follows: In section 2 we briefly present a brief overview of all of the models. Next, in section 3 we present our benchmark PC-VCSEL structure and show the simulation results obtained by each method. We also show an in-depth analysis of the case where different approaches show qualitatively different results and propose possible explanations of these discrepancies. Finally we conclude about the limits of the presented models and future challenges of precise and efficient modeling of PC-VCSELs.
2. Computational Methods
In our research we compare the results obtained with four different methods. Two of them are wholly vectorial ones and are capable of considering all the properties of light in complex devices. The third scalar method does not consider the polarization properties of the light, but has been confirmed to provide comparable results to vectorial solutions. The forth one is a hybrid vectorial-scalar model, which can distinguish different light polarizations, although it separates vectorial Helmholtz equation into two parts solved individually for different coordinates in a manner similar to the effective-index method [9]. The particulars of each of the methods are skipped from this paper. However, brief descriptions are presented below to familiarize the reader with the general ideas behind the simulations.
2.1. Coupled Mode Model (CMM)
The coupled mode model (CMM) for VCSEL simulations was suggested by P. Debernardi and his colleagues in 2001 [10]. In this model, the electromagnetic (EM) field is described by a modal expansion. As schematically illustrated in Fig. 2, a VCSEL structure is modeled to be filled with a uniform reference medium. Thus, the basis modes are free space modes of the reference medium and are common throughout all layers of the entire VCSEL structure. The EM field propagation is described by coupling of basis modes originating from perturbations of the refractive index compared to the reference medium; each different layer results in different coupling coefficient matrix. Since the basis modes are common throughout the entire structure, no mode matching conversion is required at layer interfaces, which involves heavy matrix inversion calculations especially for complex geometries. Employing free space modes as basis mode also removes transverse boundary conditions such as metallic boundary conditions with absorbing material coatings. Two longitudinal boundary conditions can be imposed typically at topmost and bottommost interfaces of the VCSEL structure, as shown in Fig. 2(b). This gives an eigenvalue equation. By solving the equation and requiring a self consistency after a round trip propagation of the field over the VCSEL structure, the lasing wavelength and threshold material gain of a mode are obtained. The field profile of the mode can be constructed from the eigenvector, which is a list of basis mode expansion coefficients.
Fig. 2. (a) A schematic profile of a VCSEL structure. (b) A modeling of the structure of (a) in CMM.
Employing modal expansion rather than spatial discretization and the absence of mode matching conversion at layer interfaces leads to very high computational efficiency in terms of simulation time and memory requirement. Finding a mode of the photonic crystal VCSEL investigated in this paper took 15–30 minutes with 100M-300MB memory allocation for an Intel 3.00 GHz E6850 processor.
2.2. Finite Element Method (FEM)
Finite Element Method is used to solve either vectorial or scalar Helmholtz equation in which the unknowns are the electric field profile E(r) and complex angular frequency (ω). The imaginary part of the angular frequency gives the rate at which the intensity of the laser mode attenuates in the lack of pumping. The threshold gain can be estimated by dividing the cold-cavity modal loss with the total confinement factor for the quantum wells. However, if the built-in mode confinement is very low and significant gain-guiding occurs near the threshold condition, the cold-cavity model is not appropriate. An estimated threshold gain is then incorporated into the complex refractive index of the active region, and the obtained modal loss is used as a correction to the assumed threshold gain.
To obtain a numerical solution one selects an appropriate computational domain, which must be covered by non-reflecting boundary condition. This task is performed by perfectly matched layers (PML), which involve a complex diagonal tensor [11]. If the unknown electric field is approximated as a linear combination of appropriately selected vectorial basis functions, E(r) = ∑i ciEi(r), an algebraic eigenproblem can be derived from the equivalent variational functional:
where c denotes the column vector of the unknown coefficients, A and B matrices can be found elsewhere [12].
The computational window is filled by prism elements. A triangular mesh is first created, which respects all lateral contours present in the PC-VCSEL. To this end, the electrical aperture, etched hole patterns, metallic contacts, and also the artificial PML region is projected into the same transverse plane. The triangularization of this domain results in a two-dimensional mesh (see Fig. 3), which is then extruded to the three-dimensional structure. Note that only a quarter cross-section is meshed, because symmetry boundary conditions can be enforced at the side surfaces [12].
Fig. 3. Sample triangular lateral mesh for a single-defect PC-VCSEL. In order to exploit symmetry conditions, only a quarter cross-section is meshed for vectorial, and a 30-degree section for scalar solutions. Etched regions are shown with blue colors, and can be matched to various hole diameters. Finer mesh is created around the defect to focus on the most important part of the laser.
The vectorial Helmholtz equation can be rightfully approximated with its scalar version, if there are no sharp material discontinuities in the laser resonator. Although this condition is not fulfilled in PC-VCSELs, the scalar solution predicts correctly the tendencies for the mode profiles, modal wavelengths and optical losses at much lower computational cost. It is enough to consider a 30-degree section to calculate the fundamental mode. The computational window is also filled with prism elements, and the scalar unknowns are associated with vertices (or the mid-points of vertical edges) of the prisms. The finite volume discretization method is described in details elsewhere [13], and leads to a similar problem to Eq. (1).
Eight layers of prisms are assigned to each quarter-wavelength layer, and selecting 40 divisions along the radius of the electrical aperture results in about 2.5 million unknowns. The calculation of a single laser mode took about half a day on 2.4 GHz Opteron processor, memory cost was around 15–20 GB. The vectorial method has been applied only for the simulation of few PC-VCSELs to validate the scalar approximation, although in much coarser resolution due to limited memory resources.
2.3. Plane-Wave Admittance Transfer (PWAT)
Plane-wave admittance transfer method (which is described in details in Refs. [14] and [15]) is capable of determining eigenmodes of any planar photonic structure. It solves Maxwell equations in a frequency domain by using a plane-wave expansion within each layer and an analytical solution in the perpendicular direction. Such approach is possible due to transformation of the EM field to the diagonal coordinates, in which matrix Helmholtz equations of the form
are transformed into a set of decoupled equations
where Γ2 is diagonal. The fields are matched using the admittance transfer technique and the eigenmode if found by searching for a complex frequency that yields a non-zero EM field distribution. In order to obtain a laser threshold gain, the imaginary part of the quantum-well refractive index is modified in each step of the root-finding algorithm, so that the imaginary part of the eigenfrequency is equal to zero.
The lateral field distribution is approximated by the truncated plane-wave expansion, which strongly reduces the dimensions of the used matrices as compared to the spatial finite-difference probing. The method is extremely efficient for modeling periodic structure due to the natural periodic boundary conditions. For analysis of single devices it is necessary to include PMLs as absorbing boundaries. This is straightforward as the method directly considers anisotropic materials with complex permittivity and permeability.
The exact determination of each mode is performed with the Broyden root-finding algorithm. An average computation of a single mode requires around 1.5GB of memory and takes approximately 4 hours on a single Intel Xeon 2.5 GHz processor.
2.4. Plane-Wave Effective Index (PWEI)
The plane-wave effective index method is a simplified version of the plane-wave admittance transfer method. It considers a planar structure and expands the lateral field distribution in plane-wave basis. However, the electric field is separated into two components, which are determined separately for directions parallel and perpendicular to the device layers
The difference of this method compared to the popular effective-index technique is the fact that the EM field at all times is considered as a vectorial one, hence the polarization information is not lost [16]. The eigenmodes are determined as complex eigenvalues of a two-dimensional problem, which gives a significant performance boost with respect to computation time and memory requirements, while still enabling computation of modal losses and threshold gain. However, due to the application of Eq. (4) some properties of light in more complex structures are not considered properly, especially in cases where the mode is not strongly confined [16].
The main advantage of this approach is its efficiency. The computations of a a single mode took less than 1 minutes on 2.5 GHz Intel Xeon with memory consumption of 400MB. The PWAT and PWEI methods are implemented in an open-source software PSlab (http://pslab.sourceforge.net).
3. Results
The structure analyzed in our benchmark is a gallium arsenide PC-VCSEL presented in Fig. 1. It consists of single-wavelength long GaAs cavity, sandwiched by 24/29 pairs of top/bottom GaAs/AlGaAs distributed Bragg reflectors (DBRs). The optical gain is provided by three layers of 8 nm-thick InGaAs quantum wells separated by 5 nm GaAs barriers and located together in the anti-node position. The diameter of the optical gain region is set to be the same as the photonic crystal pitch Λ. Details for layer thickness and refractive index values can be found in Tab. 1.
Table 1. Details of the layer structure of the analyzed VCSEL.
The photonic crystal is etched in the structure and consists of three rings of a triangular lattice with one missing hole in the center. In the analysis four variants are considered: (i) the holes are etched through the whole laser structure, (ii) the holes are present in the cavity only, while DBRs are uniform, (iii) the holes are etched in the whole top DBR, and (iv) only 10 pairs of the top DBR are etched. The PC is the only factor providing light confinement in the VCSEL. Without it, the resonance wavelength would be 980 nm, with the vertical cavity Q-factor around 20 000.
In the present benchmark, we focus on the fundamental mode of the PC-VCSEL. Due to the Gaussian-like shape, this mode is the most preferred one in most applications. Furthermore VCSELs with incorporated photonic crystal are known for its high mode selectivity [3], which makes it feasible to tune them into a single-mode source with high power emission.
Fig. 4. Resonant wavelengths of PC-VCSELs with photonic crystals etched in different parts of the structure as a function of the hole diameters.
Our results are presented in Figs. 4 to 7. Figs. 4 and 5 depict computed resonant wavelengths (λ) and threshold gains (g) as functions of the hole diameters for the photonic-crystal pitch Λ equal to 4.0 µm. The pitch of 4.0 µm is a typical value and has been kept fixed for the above mentioned cases (i) to (iv). In Figs. 6 and 7 we show the results determined for different sizes of the photonic-crystal lattice constant for cases (i) and (iii). From these graphs one can observe that all of the presented numerical methods show similar tendencies in the simplest cases, however, there are some quantitative discrepancies. Furthermore there exist some qualitative differences in some of the plots.
Fig. 5. Threshold gains of PC-VCSELs with photonic crystals etched in different parts of the structure as a function of the hole diameters.
First of all the FEM method showed a systematic shift in the wavelength (not present in the figures), which was the effect of the finite axial resolution. This shift was similar for all configurations and it did not affect the more important wavelength spacing. Hence, the necessary correction was determined by comparing the 1D simulations done with the usual λ/32 and much finer λ/10000 resolutions. The obtained value of 1.367 nm was subtracted from all the FEM wavelengths. After such operation, all the analyzed methods predict the same resonant wavelength with good agreement. The only exception is PWEI for the most shallow holes (Fig. 4d), in which case a visible drop of λ for increasing diameter of holes is observed. This behavior can be explained by the fact that this drop is similar to the one predicted by all methods for deeper holes (Figs. 4a–c). The shallow holes do not differ qualitatively from deep ones for the effective index approach (i. e. the only difference is in the effective index in the hole-region projected into the x-y plane).
Unsurprisingly, stronger differences are present in case of the threshold gain data. All the methods show good agreement only in case of the Λ = 4µm wholly etched structure. All other situations require more detailed analysis. Fig. 5c shows the threshold gains for holes etched only in the top DBR. Due to their complicated arrangement, large scattering occurs at their edges, which is reflected in an increase in the threshold gain for d/Λ > 0.4. This effect is not predicted by PWEI, which—as a variation of the effective index method—is incapable of properly considering scattering losses. Similarly, it even stronger underestimates the threshold gain for shallow holes (Fig. 4d), where scattering is much larger due to the weaker mode confinement. Furthermore, this weak confinement forces an à priori introduction of an estimated threshold gain in the quantum wells for FEM computations, as described in section 2.2. The same procedure is utilized by FEM also in other configurations, like for d/Λ = 0.1.
Fig. 6. Resonant wavelengths of PC-VCSELs with photonic crystals etched in different parts of the structure as a function of the photonic crystal pitch.
Fig. 7. Threshold gains of PC-VCSELs with photonic crystals etched in different parts of the structure as a function of the photonic crystal pitch.
If the holes are etched only in the cavity, the mode is well confined as in the case of the wholly-etched structure. However, CMM and FEM show a significant increase of the threshold gain for holes with diameters larger than 0.3Λ. This effect is not represented by both PWEI and PWAM, although the latter one shows better agreement with CMM and FEM for increased number of planewaves: the threshold gain for d/Λ = 0.8 raises to 570 cm−1, which is much closer value to the two other methods. However, in such a case the requirement for computer resources and simulation times increases so strongly that we decided not to perform additional simulations with even higher numbers of plane waves. In order to investigate the origin of the predicted threshold gain by the different models, we analyzed the vertical (i. e. in a plane perpendicular to the VCSEL layers) cross-sections of the electric field profiles computed with CMM, FEM and PWAT. They are presented in Fig. 8. As one can see, CMM and FEM predict diagonally directed fields originating near the air holes in the cavity, which correspond to the diffraction scattering. However, in case of the PWAT with the default number of plane-waves, this effect is much weaker, especially when compared with the veritcally emitted flux.
Fig. 8. Electric field intensity profiles in the vertical cross-section of the PC-VCSEL, computed for Λ = 4µm, d/Λ = 0.5 with holes etched in the cavity only. (a) xyz profile in dB scale obtained using the CMM, (b) xz profile in log scale obtained using the FEM, and (c) xz profile in dB scale obtained using PWAT.
We can observe a reasonable agreement between all the methods for the wholly etched structure (Fig. 7a), when we vary the photonic-crystal lattice constant for a fixed d/Λ ratio set to 0.5; their mutual gain difference does not exceed 17%. It is noteworthy that for Λ larger than 6µm the CMM is unable to find accurate solutions. It is because the coupling coefficient matrix becomes close to singular or badly spaced when the total air hole volume is very large. This singular or badly-spaced matrix is difficult to be accurately inversed or solved for eigenvalues, using Matlab. However, we may say that for typical PC-VCSEL designs, the total air hole volume is much smaller. Here, the CMM can accurately anticipate the threshold gain. Fig. 7b shows the threshold gains for the holes etched only in the top DBR. Contrary to the wholly etched case in Fig. 7a, the PWEI shows considerable discrepancy from the other methods, for small pitches. This is caused by poor mode confinement and strong scattering for small Λ. Therefore, the correct threshold gain can not be predicted by effective index approaches in these cases, including PWEI.
4. Conclusions
In this article we have presented the comparison of the computed resonant wavelength and threshold gain of a PC-VCSEL computed with four different numerical models. All the methods showed a very good agreement in the wavelength, apart of the PWEI for the case of shallow holes. The computed threshold gain was also predicted with reasonable agreement except some exceptional cases: very large hole volume for CMM, small hole diameter for FEM and holes etched only in the cavity for PWAT and PWEI (which both shared the plane-wave expansion for the representation of the lateral distribution of the field). Furthermore, the PWEI yielded different results in some other cases, which clearly demonstrated the limitations of the effective-index approximations. The main reason of the shown discrepancies was the three-dimensional complexity of the structure and strong scattering of the light at the edges of the holes. Thus, it was difficult to reach similar level of agreement between methods as achieved in the benchmark performed in the past for a classical oxide-confined VCSEL with axial symmetry [17]. None of the methods is perfect and one must be aware of their limitations when performing numerical simulations of such complex structures as PC-VCSELs. Although in all the cases, the differences between CMM, FEM and PWAT could be reduced by using finer mesh or larger number of basis functions, such procedure would significantly increase computation time and required resources.
In conclusion, our benchmark has clarified issues that can arise in optical modeling of static properties of PC-VCSELs. In case of dynamic simulations, which nowadays gain more and more interests, the situation will not be less complicated. Additional factors would play an important role like for example non-linear effects, which will limit the precision of the numerical simulations. However, if one is aware of the computational limitations of the different methods, it is possible to choose the one which best fits the structure under investigation.
Appendix
The numerical values of our results are presented in tables for the readers reference. Tab. 2 shows the data computed for a constant pitch Λ = 4µm and varied hole radius and Tab. 3 refers to the computations done for varied pitch and the constant d/Λ = 0.5 ratio.
Table 2. The results summary computed for Λ = 4.0µm.
Table 3. The results summary computed for d/Λ = 0.5.
Acknowledgments
We would like to acknowledge the European COST Action MP0702, which made the presented investigations of the benchmark structure possible. COST states for European Cooperation in Science and Technology. M.D. acknowledges the Polish Ministry of Science for a grant no. 505/N-COST/2009/10.
References and links
1. D. S. Song, S. H. Kim, H. G. Park, c. K. Kim, and Y. H. Lee, “Single-fundamental-mode photonic-crystal vertical-cavity surface-emitting lasers,” Appl. Phys. Lett. 80, 3901–3903 (2002). [CrossRef]
2. A. J. Danner, T. S. Kim, and K. D. Choquette, “Single fundamental mode photonic crystal vertical cavity laser with improved output power,” Electron. Lett. 41, 325–326 (2005). [CrossRef]
3. T. Czyszanowski, M. Dems, and K. Panajotov, “Single mode condition and modes discrimination in photonic-crystal 1.3 µm AlInGaAs/InP VCSEL,” Opt. Express 15, 5604–5609 (2007). [CrossRef] [PubMed]
4. T. S. Kim, A. J. Danner, D. M. Grasso, E. W. Young, and K. D. Choquette, “Single fundamental mode photonic-crystal vertical cavity surface emitting laser with 9 GHz bandwidth,” Electron. Lett. 40, 1340–1341 (2004). [CrossRef]
5. S. Bischoff, F. Romstad, M. Juhl, M. H. Madsen, J. Hanberg, and D Birkedal, “2.5 Gbit/s modulation of 1300 nm single-mode photonic crystal VCSELs,” in “Optical Fiber Communication Conference and Exposition and The National Fiber Optic Engineers Conference, Technical Digest (CD),” (2006), OFA6.
6. F. Romstad, S. Bischoff, M. Juhl, S. Jacobsen, and D. Birkedal, “Photonic crystals for long-wavelength single-mode VCSELs,” Proc. SPIE 6908, 69080C (2008). [CrossRef]
7. D. S. Song, Y. J. Lee, H. W. Choi, and Y. H. Lee, “Polarization-controlled, single-transverse-mode, photonic-crystal, vertical-cavity, surface-emitting lasers,” Appl. Phys. Lett. 82, 3182–3184 (2003). [CrossRef]
8. M. Dems, T. Czyszanowski, and K. Panajotov, “Highly birefringent and dichroic photonic-crystal VCSEL design,” Opt. Commun 281, 3149–3152 (2008). [CrossRef]
9. G. R. Hadley, “Effective index model for vertical-cavity surface-emitting lasers,” Opt. Lett. 20, 1483–1485 (1995). [CrossRef] [PubMed]
10. G. P. Bava, P. Debernardi, and L. Fratta, “Three-dimensional model for vectorial fields in vertical-cavity surface-emitting lasers,” Phys Rev. A 63, 23816 (2001). [CrossRef]
11. S. D. Gedney, “An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices,” IEEE T. Antenn. Propag. 44, 1630–1639 (1996). [CrossRef]
12. P. Nyakas, “Full-vectorial three-dimensional finite element optical simulation of vertical-cavity surface-emitting lasers,” IEEE J. Lightwave Techn. 25, 2427–2434 (2007). [CrossRef]
13. P. Nyakas, G. Varga, Z. Puskás, N. Hashizume, T. Kárpáti, T. Veszprémi, and G. Zsombok, “Self-consistent real three-dimensional simulation of vertical-cavity surface-emitting lasers,” J. Opt. Soc. Am. B 23, 1761–1769 (2006). [CrossRef]
14. M. Dems, R. Kotynski, and K. Panajotov, “Plane-wave admittance method — a novel approach for determining the electromagnetic modes in photonic structures,” Opt. Express 13, 3196–3207 (2005). [CrossRef] [PubMed]
15. M. Dems, “Plane-wave admittance method and its applications to modeling semiconductor lasers and planar photonic-crystal structures,” Ph.D. thesis, Technical University of Lodz (2007).
16. M. Dems and K. Panajotov, “Modeling of single- and multimode photonic-crystal planar waveguides with plane-wave admittance method,” Appl. Phys. B 89, 19–23 (2007). [CrossRef]
17. P. Bienstman, R. Baets, J. Vukusic, A. Larsson, M. J. Noble, M. Brunner, K. Gulden, P. Debernardi, L. Fratta, G. P. Bava, W. Wenzel, K. Klein, C. Conradi, R. Pregla, S. A. Riyopoulos, J.-F. P. Seurin, and S. L. Chuang, “Comparison of optical VCSEL models on the simulation of oxide-confined devices,” IEEE J. Quantum Electron. 37, 1618–1631 (2001). [CrossRef]
