Design of nanowire-based vertical-cavity surface-emitting lasers Master’s thesis in Master Programme Applied Physics MARCUS BENGTHS Photonics Laboratory, Department of Microtechnolgy and Nanoscience CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2018 Master’s thesis 2018 Design of nanowire-based vertical-cavity surface-emitting lasers MARCUS BENGTHS Photonics Laboratory Department of Microtechnolgy and Nanoscience – MC2 Chalmers University of Technology Gothenburg, Sweden 2018 Design of nanowire-based vertical-cavity surface-emitting lasers MARCUS BENGTHS © MARCUS BENGTHS, 2018. Supervisors: Assoc. Prof. Åsa Haglund Assoc. Prof. Johan Gustavsson PhD student Filip Hjort Assoc. Prof. Jörgen Bengtsson Examiner: Assoc. Prof. Åsa Haglund Master’s Thesis 2018 Photonics Laboratory Department of Microtechnology and Nanoscience – MC2 Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Norm of electrical field inside platelet-based GaN VCSEL with concave top DBR and platelet radius of 1300 nm. See Section 4.1 for further information. Typeset in LATEX Printed by Chalmers Tryckeri – Reproservice Gothenburg, Sweden 2018 iv Design of nanowire-based vertical-cavity surface-emitting lasers MARCUS BENGTHS Photonics Laboratory, Department of Microechnology and Nanoscience Chalmers University of Technology Abstract The optical losses have been studied in a novel nanowire-based GaN vertical-cavity surface-emitting laser as a function of different design parameters. The simulations were based upon the finite element frequency domain method and were implemented in COMSOL Multiphysics® simulation software. The results show that the thresh- old material gain depends most strongly on the platelet radius. Lasers with a small platelet radii are shown to suffer from large and inevitable diffraction losses. Hence, in order to reach decent threshold material gain and achieve lasing with the cur- rent geometric design, the platelet radius must be increased to at least 1200 nm to 1300 nm. Additional diffraction losses derive from the pitched top of the platelet, but these can be suppressed by depositing a concave top distributed Bragg reflector (DBR) with focusing abilities. The threshold material gain can be further reduced by some extent by adding additional mirror pairs in the top DBR. Furthermore, the thickness and materials of the cladding layer also strongly affects the threshold material gain. However, the influence of the cladding layer is clearly different be- tween having a flat top DBR or a concave top DBR. This complex influence of the cladding layer requires further investigation. The method presented has been proven robust, and comparable results to a beam propagation method (BPM) validates its accuracy. Being far less computationally heavy than BPM, it will be an important tool in the future design of novel lasers. Keywords: VCSEL, GaN, COMSOL Multiphysics®, FEFD v Acknowledgements First, I would like to express my deep gratitude to my main supervisor and examiner Assoc. Prof. Åsa Haglund for her enthusiastic encouragement and rapid feedback. I would like to thank Assoc. Prof. Johan Gustavsson and Assoc. Prof. Jörgen Bengtsson for their guidance and fruitful discussions that enabled the development of the threshold material gain extraction method. I am also grateful to Filip Hjort for the introduction to COMSOL Multiphysics® that gave me a solid start with my thesis. I would like to extend my thanks to the great employees at the Photonics Laboratory that have made the lunch discussions most joyful. Finally, I would like to thank my family for always supporting me and my friends for making my time at the university to memorable time of my life. Marcus Bengths, Gothenburg, June 2018 vii Contents List of Tables xi List of Abbreviations xiii 1 Introduction 1 1.1 Aim . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Theory 3 2.1 Working principles of a semiconductor laser . . . . . . . . . . . . . . 3 2.1.1 Conventional VCSEL Structure . . . . . . . . . . . . . . . . . 6 2.1.2 Nanowire-based GaN VCSEL Structure . . . . . . . . . . . . . 7 2.2 Maxwell’s Equations and the wave equation . . . . . . . . . . . . . . 8 2.2.1 Frequency domain representation . . . . . . . . . . . . . . . . 9 2.2.2 Physical assumptions . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Optical Modes in VCSELs . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.1 Longitudinal Modes . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.2 Transversal Modes . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.3 Cylindrical geometry approximation . . . . . . . . . . . . . . . 11 2.4 The Finite Element Frequency Domain Method . . . . . . . . . . . . 12 3 Methods 15 3.1 COMSOL Multiphysics® . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.1.1 Geometries of interest . . . . . . . . . . . . . . . . . . . . . . 15 3.1.2 Implementing Physics . . . . . . . . . . . . . . . . . . . . . . 17 3.1.2.1 Boundary Conditions . . . . . . . . . . . . . . . . . . 18 3.1.2.2 Perfectly Matched Layers (PML) . . . . . . . . . . . 19 3.1.3 Meshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.1.4 Solver Configuration . . . . . . . . . . . . . . . . . . . . . . . 21 3.2 LivelinkTM for MATLAB® . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3 The Main Program . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3.1 Gain Sweep Algorithm . . . . . . . . . . . . . . . . . . . . . . 24 3.3.1.1 Threshold correction . . . . . . . . . . . . . . . . . . 26 3.4 Post-Processing Algorithms . . . . . . . . . . . . . . . . . . . . . . . 26 3.4.1 List-shift detection . . . . . . . . . . . . . . . . . . . . . . . . 27 3.4.2 List-shift processing . . . . . . . . . . . . . . . . . . . . . . . . 28 ix Contents 4 Results 31 4.1 Platelet Radius . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2 Cladding thickness and refractive index . . . . . . . . . . . . . . . . . 35 4.3 Concave DBR structure . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.3.1 Cladding thickness and refractive index . . . . . . . . . . . . . 39 4.4 Number of DBR mirror pairs . . . . . . . . . . . . . . . . . . . . . . 40 4.5 Active region alignment . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.6 Absorption effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.7 Reliability and error estimates . . . . . . . . . . . . . . . . . . . . . . 43 4.7.1 Simulation Box Radius . . . . . . . . . . . . . . . . . . . . . . 44 4.7.2 PML Thickness . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.7.3 Maximum Element Size . . . . . . . . . . . . . . . . . . . . . 46 5 Conclusions 47 Bibliography 49 A Derivation of transversal modes of a step-index optical fiber I x List of Tables 3.1 List of Geometrical Parameters . . . . . . . . . . . . . . . . . . . . . 16 3.2 List of Vertical Parameters . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3 List of Radial Parameters . . . . . . . . . . . . . . . . . . . . . . . . 17 3.4 List of Material Parameters . . . . . . . . . . . . . . . . . . . . . . . 18 xi List of Tables xii List of Abbreviations Abbreviation Meaning AlxGa1−xN Aluminum-gallium-nitride AlN Aluminum-nitride CAD Computer-aided design DBR Distributed Bragg reflector EBL Electron-blocking layer EEL Edge-emitting laser FDTD Finite difference time domain FEFD Finite element frequency domain FEM Finite element method GaAs Gallium-arsenide GaN Gallium-nitride HEmp Hybrid mode InxGa1−xN Indium-gallium-nitride InN Indium-nitride ITO Indium tin oxide LED Light-emitting diode MOCVD Metal-Oxide chemical vapour deposition MQW Multiple quantum well n-GaN n-doped gallium-nitride NW Nanowire p-GaN p-doped gallium-nitride PML Perfectly matched layer QW Quantum well QWb Quantum barrier section in MQW structure QWw Quantum well section in MQW structure SAG Selective area growth SBC Scattering boundary conditions Si3N4 Silicon nitride SiO2 Silicon dioxide TE0p Transversal electric mode TM0p Transversal magnetic mode UV Ultraviolet VCSEL Vertical-cavity surface-emitting laser VCS Vertical cross-section xiii List of Tables Parameter Meaning a Side length of platelet hexagon A Unknown prefactor αabs Absorption losses αd Diffraction losses αGaN Absorption coefficient for GaN αITO Absorption coefficient for ITO αm Mirror losses B Unknown prefactor ~B Magnetic field density β Propagation constant c0 Speed of light in vacuum C Unknown prefactor dair Air region thickness dDBR,high High index DBR layer thickness dDBR,low Low index DBR layer thickness dEBL EBL thickness dITO ITO (indium tin-oxide) thickness dmask Si3N4-mask thickness dnGaN n-GaN length dnGaNp n-GaN length (pitched only) dnGaNv n-GaN length (vertical only) dpl Platelet height dQWb Quantum well barrier thickness dQWw Quantum well well thickness dpGaN1 p-GaN length (below EBL) dpGaN2 p-GaN length (above EBL) dsapphire Sapphire substrate thickness dspacing Spacing layer thickness dsubstrate GaN substrate thickness D Unknown prefactor ~D Displacement field δ Damping of the solution/ damping factor δ0 y-intercept in damping factor and applied gain relation δactive Active region displacement from default δi Damping of the solution from gain step i δli Damping of the solution from gain step i, solution index l. δ̃ Actual damping of the solution δ̃0 Actual y-intercept in damping factor and applied gain relation δ̄ Mean value of damping of solution for one gain sweep |∆λ| Longitudinal mode separation ∆ω Thickness of PML in COMSOL definition EC Conduction band minimum Eg Bandgap Ei Electric field expansion coefficient in test function basis Eφ Azimuthal component of electric field Er Radial component of electric field xiv List of Tables Parameter Meaning EV Valence band maximum Ez Vertical component of electric field ~E Electric field strength ε0 Vacuum permittivity [εr] Relative permittivity tensor η Arbitrary composant index fν Load vector fref Reference eigenfrequency (frequency guess) f(ξ) Complex PML stretching function f ′(ξ) Real part of PML stretching function f ′′(ξ) Imaginary part of PML stretching function F (x) PML analytical continuation function φ Azimuthal angle coordinate φ̂ Azimuthal base vector Φ(φ) Azimuthal separation function Φ′′(φ) Second derivative of azimuthal separation function g Applied optical material gain gi Applied gain at gain step i gth Threshold material gain g̃ Actual applied gain g̃th Actual threshold material gain ḡ Mean value of applied gain for one gain sweep γ Cladding shifted propagation constant Γactive Optical confinement factor h Maximum element size Hφ Azimuthal component of magnetic field Hr Radial component of magnetic field Hz Vertical component of magnetic field ~H Magnetic field intensity i Gain sweep step index j Imaginary unit Jm m-th order Bessel function of first kind J ′m First derivative of m-th order Bessel function of first kind ~J Current density k Slope between damping factor and applied gain k2 0 Square magnitude of vacuum wavevector kli Local slope between damping factor and applied gain kx 1D wavevector k̃ Actual slope between damping factor and applied gain Km m-th order modified Bessel function of second kind K ′m First derivative of m-th order modified Bessel function of second kind Kµν Stiffness matrix κ Extinction coefficient l Eigenvalue index (Solution index) L Cavity length Leff Effective cavity length xv List of Tables Parameter Meaning λ Wavelength λ0 Vacuum wavelength λDBR Center wavelength of DBR reflectivity peak λref Wavelength guess λt PML typical wavelength λ̃ Angular eigenfrequency m Azimuthal wavenumber M Denominator in maximum element size expression Mµν Mass matrix µ Degree of freedom index µ0 Vacuum permeability [µr] Relative permeability tensor n Refractive index nair Refractive index of air ncore Refractive index of core nclad Refractive index of cladding nDBR,high Refractive index of high index material in DBR nDBR,low Refractive index of low index material in DBR neff Effective refractive index nembedding Refractive index of Embedding Material nEBL Refractive index of EBL nGaN Refractive index of GaN nHfO2 Refractive index of HfO2 in DBR nITO Refractive index of ITO nlocal Local refractive index nmask Refractive index of Mask nsapphire Refractive index of Sapphire Substrate nSiO2 Refractive index of SiO2 in DBR nQWb Refractive index of QW barrier nQWw Refractive index of QW well N Number of anti-nodes inside the cavity Ndof Total number of degrees of freedom Neigs Number of requested eigenvalues NGS Number of gain sweep steps ~Nµ Vector test function ν Degree of freedom index p Mode solution index P (r) Radial separation function P ′(r) First derivative of radial separation function P ′′(r) Second derivative of radial separation function Π PML curvature factor Q Sum of the square of the residual r Radial coordinate r̂ Radial base vector ~r Space vector xvi List of Tables Parameter Meaning RClad Cladding layer thickness Rembedding Embedding layer width Rmask Mask hole radius Rpl Platelet radius RplTop Platelet radius at top end RPML Thickness of PML Rsimbox Simulation Box Size s PML scaling factor [σ] Conductivity tensor 2σ Double standard deviation σl Standard deviation of local slopes for solution index l t Time θ Tilt angle on platelet top section V Volume of interest ∂V Boundary of volume of interest x x-coordinate in Cartesian coordinate system x0 x-coordinate for physical domain - PML interface xedge x-coordinate for edge of simulation box ξ Normalized PML coordinate z Vertical coordinate ẑ Vertical base vector Z(z) Vertical separation function Z ′′(z) Second derivative of vertical separation function ζ Core shifted propagation constant ω Angular frequency xvii List of Tables xviii 1 Introduction There are two main types of semiconductor lasers. Those are surface-emitters and edge-emitters. The vertical-cavity surface-emitting laser (VCSEL) is a type of semi- conductor laser architecture where the laser beam is emitted vertically in relation to the substrate, in contrast to the commonly known edge-emitting laser (EEL) which emits light horizontally. The VCSEL has several advantages compared to EEL. That includes low threshold currents, circular, low-divergent output beam, and low-cost manufacturing [1]. In addition, since the VCSEL emits light vertically it is possible to directly make use of the 2D-array structure. Today, VCSELs are found in a diverse field of applications such as computer mice, laser printers, sensors, and illumination, to name a few. Most of the com- mercial VCSELs are based on GaAs and emitt light in the red to infrared regime (∼ 800 nm), thus there is a lack of commercially viable VCSELs operating in the short wavelength part of the optical spectrum [1]. Hence, for the past decade, much attention has been admitted to III-nitride materials, i.e. (Al, In, Ga)N, which can be used for light-emitting devices operating in the blue-ultraviolet regime (∼ 400 nm). Such device has many possible applications, e.g. high-resolution printers and sen- sors, visible light communication, high-density optical data storage, and to replace blue in white light emitting diodes (LEDs) thus providing more power-efficient illu- mination at high drive current densities in for example displays. It can also be used in various medical fields, such as disinfection and caries-detection within dentistry, optogenetics, and cancer diagnosis[1]. However, although much research has been done during the past decade, a power- and cost efficient blue-UV GaN-lasers is yet to be developed, thus it implies that this is a challenging task. A fundamental issue when composing III-nitride material heterostructures, e.g. distributed Bragg Reflectors (DBRs) and quantum well structures, is that it is dif- ficult to find materials with sufficient refractive index contrast that also are lattice- matched. A lattice mismatch induces stresses in the system, which can result in high dislocation densities, which are defects that strongly limits the performance of the laser [1]. A promising alternative to planar GaN lasers is that of nanowire-based lasers. Nanowires are due to their geometry free to relax in plane. Therefore lattice mismatch is no longer of concern, and it is possible to grow nanowire cavities that are both strain- and dislocation-free [2]. This Master Thesis is a part of an ongoing research project on a new type of nanowire-based GaN VCSEL. The project is a collaboration between Chalmers University of Technology and Lund University. The nanowires are grown by selective area growth (SAG) by Metal-Oxide Chemical Vapour Deposition (MOCVD), by Lars Samuelson’s group at Lund University. They use a novel growth method that 1 1. Introduction makes it possible to grow large dislocation-free GaN nanowires with diameters up to 1000 nm, henceforth referred to as ”platelets”. Åsa Haglund’s group at Chalmers University of Technology has developed an electrochemical lift-off technique by which it is possible to achieve smooth surfaces, which allows deposition of both top and bottom dielectric DBR onto the device. By combining these two novel techniques, we hope to be able to realize an electrically pumped platelet-based GaN VCSEL that outperforms today’s blue-UV emitters in both power efficiency, size and cost- efficiency. 1.1 Aim This master thesis focus on optical simulations of the laser cavity, for the purpose of understanding the loss-mechanisms of this particular type of laser structure, and to propose good laser designs. The optical simulations will be be based on finite- element frequency domain techniques using COMSOL Multiphysics® simulation software. Optical mode field distributions, transverse optical confinement, resonator losses, and threshold gain will be determined for different geometrical and material compositional variations of the structure. Variations include design of cladding layers, planarity and width of distributed Bragg reflectors (DBRs), and platelet geometries. Based on these simulation results, laser designs that improve laser performance will be proposed. 1.2 Scope This thesis will mainly focus on optical simulations of the laser structure. Electrical properties such as charge carrier distributions and contact resistance won’t be cov- ered, nor will be the inclusion of thermal effects on optical and electrical properties, but those are kept in mind when proposing new structural designs. Fabrication and realization of the devices is beyond the scope of this master thesis, but the new designs will be discussed with experts in these fields. 2 2 Theory This chapter serves as a theoretical walkthrough for concepts relating to this project. The first part is an introduction to the fundamental working principles of a semicon- ductor laser, first in general, and later on applied to vertical cavity surface emitting lasers (VCSEL). An overview of the platelet-based VCSEL is presented as well. The second part focuses on the derivation of the electromagnetic wave equation and highlights some physical assumptions made in the modelling. The third part de- scribes the concept of longitudinal and transversal optical modes, and ends with an important comparison between cylindrical and hexagonal platelets. The last section is a brief introduction to the finite element method and the finite element frequency domain techniques that is the main numerical method used in the simulations pre- sented in this thesis report. 2.1 Working principles of a semiconductor laser Laser is an acronym that stands for Light Amplification by Stimulated Emission of Radiation. It is type of light that is characterized as closely monochromatic and highly coherent. A semiconductor laser is a laser device made of semiconductor materials. Let’s consider a two band material system with a valence band maximum EV and a conduction band minimum EC , where the band gap is Eg = EC −EV . As light interacts with such system, there are three processes that can occur (See Figure 2.1). In the presence of a photon with sufficient energy, the photon may be absorbed by the material stimulating an excitation of an electron from the valence band to the conduction band. This process is called electromagnetic absorption (I). The opposite process is that of spontaneous emission (IIa), which is when an excited electron transition from the conduction band to the valence band, where as the excessive energy is emitted in form of a photon. The third type of process is less intuitive, but serves an essential role for the existence of laser phenomena. As predicted by Einstein in 1917[3], there is a second type of emission that only occurs in presence of a photon. Besides photon absorption, the photon may alternatively stimulate deexcitation of an electron from the conduction band to the valence band, creating an exact copy of itself that has the same phase, frequency, and direction. Consequently, in comparison to the spontaneous emission (IIa), this stimulated emission (IIb) not only increases the number of photons in the system, but the emitted photons are also pairwise coherent, which is a requirement for obtaining laser light. 3 2. Theory (I) EC EV (IIa) (IIb) Figure 2.1: Illustration of (I) electromagnetic Absorption, (IIa) Spontaneous Emis- sion, and (IIb) Stimulated Emission. The transitions are here done between two energy bands EC and EV . Normally, the rate of spontaneous emission dominates over the rate of stimulated emission. Hence, most light sources based on electroluminescence emit incoherent light[4]. In order to overcome this natural hindrance and make a coherent light source based on stimulated emission, the electronic system must reach a state of population inversion. At population inversion, the probability of finding an electron in the conduction bans is higher than finding it in the valence band. Consequently, when a photon is present, the likelihood of stimulated emission is high. Hence, the number of photons in the system increases and the active medium acts as a light amplifier rather than an attenuator. This unnatural condition can be met through either applying an external optical source (optical pumping) or by increasing the number of carriers by applying a forward bias voltage (electrical pumping)[5]. Another important requirement for making a laser is to place the light amplifying material inside an optical cavity. The optical cavity serves as an optical feedback loop that helps to create and maintain a high photon density. Since a high photon density also affects the rate of stimulated emission, the light amplifying property of the active medium is enhanced. A typical laser optical cavity is formed by putting the active medium between a pair of mirrors. In addition, one of the mirrors must be semitransparent, which allows a small portion of light to be transmitted through the end of the cavity. This transmitted light is the laser beam. Due to constructive interference with reflected beams, the electromagnetic field inside the optical cavity of a laser device form a standing wave. The requirement for forming a standing wave is that the field must repeat itself after one roundtrip in the cavity. This requirement is summarized as a phase condition exp { −j 2π λ0/neff 2L } = 1⇒ λ0 = 2Lneff N , (2.1) where λ0 is the light wavelength, neff is the effective refractive index the light expe- riences as it propagates through the cavity, L is the length of the optical cavity, and N is a positive integer, which corresponds the number of anti-nodes that fits inside the cavity. Consequently, the lasing wavelength depends on the cavity structure and on the materials within. The amplification of light as it travels through the active medium is represented by an optical material gain, g. The stronger the pumping, the higher is the rate 4 2. Theory of stimulated emission, thus the higher is the gain. Unfortunately, the gain must compete with a couple of loss mechanisms. See Figure 2.2. That includes absorption losses, mirror losses, and diffraction losses. The absorption losses are dependent on the materials inside the cavity as well as the wavelength of the propagating light. Mirrors with poor reflectivities have high mirror losses. And specifically, for a device where the operating wavelengths are comparable to the structural dimensions of the cavity, diffraction effects also become significant[6]. Laser Beam (Ic)(Ic) (Ic) (Ic) (Ib) (Ia) (II) Figure 2.2: Schematic description of a typical laser optical cavity. The green sections represent semitransparent mirrors. The mirror on the right has lower re- flectivity than the mirror on the left, hence a stronger beam is transmitted on the right hand side. The roman (I) correspond to the main types of optical losses: (Ia) Electromagnetic absorption, (Ib) Mirror losses, and (Ic) Scattering losses. (II) is stimulated emission. For a sufficiently high injection rate of carriers, the optical gain will overcome these loss mechanisms and lasing will occur. The minimum gain required for lasing is called threshold material gain, which is denoted gth. The threshold condition is met when the material gain fully compensates for the net loss. Γactivegth = αabs + αm + αd (2.2) where αabs, αm, and αd stands for absorption losses, mirror losses, and diffraction losses respectively. The gain and losses are usually defined in units of inverse cen- timeters [1/cm]. The optical confinement factor, Γactive, is a fraction that describes how large portion of the total optical power that is confined within the active region of the cavity, i.e. a measure of the overlap between the active region and the optical field distribution. Hence, in order to determine the threshold material gain it is necessary to investigate the possible electromagnetic field distributions inside the particular optical cavity. 5 2. Theory 2.1.1 Conventional VCSEL Structure There are two main types of semiconductor lasers. Those are edge emitters and surface emitters. The edge emitting lasers (EEL) emitt light paralell to the active region, and vertical cavity surface emitting lasers (VCSEL) emitt light perpendicular to the active region, see Figure 2.3. As a consequence, the length in which the light is amplified is much shorter in a VCSEL compared to an EEL. Thus, in order to compensate for that, the semitransparent mirrors must have very high reflectivities. One type of structures that can exhibit sufficiently high reflectivities for use in VCSELs are distributed Bragg reflectors (DBRs)[7]. Edge Emitting Laser Vertical-Cavity Surface-Emitting Laser Figure 2.3: Comparison between edge emitters and vertical cavity surface emit- ters. The green regions represent semitransparent mirrors, and the yellow regions indicates active regions. A general DBR consists of a periodic structure with two alternating material layers. The layer thicknesses are quarter wavelength, i.e. dDBR,high = 1 4 λ0 nDBR,high , dDBR,low = 1 4 λ0 nDBR,low , (2.3) where nDBR,high is the refractive index of the high index material and nDBR,low is the refractive index corresponding to the material with the lowest index. λ0 is the vacuum wavelength. If this condition is met, the partial reflections from each layer interfere constructively and contribute to the overall mirror reflectivity. By adding more layers, the net mirror reflectivity is increased. The reflectivity of a DBR is wavelength dependent. The reflectivity spectrum is comb-shaped with a broad high reflectivity peak centered at a wavelength λDBR. Thus in order to optimize the effect of the mirrors, λDBR must match the lasing wavelength λ0. Fortunately, the reflectivity spectrum can be tuned by simply adjusting the thickness of the DBR layers. A high refractive index contrast between the DBR layers will result in a spectrally wide reflectivity peak, and a low index contrast will narrow down the reflectivity peak, which on one hand functions as a narrow band-pass filter but on 6 2. Theory the other hand makes it more difficult to match the resonance wavelength of the cavity [8]. Futhermore, as a way to reduce the threshold material gain of a VCSEL, the active region is typically made as a multiple quantum well structure (MQW) [5]. The quantum wells have a width on the order of just a couple of nanometers. Con- sequently, this confinement causes quantization of the electronic properties of the structure from 3D into 2D, and the continous conduction and the valence bands are split into subbands with discrete energies. The quantization also results in a sharper electron density distribution at the band edge, which makes it easier to acheive population inversion. Hence, the threshold material gain is lowered [5]. 2.1.2 Nanowire-based GaN VCSEL Structure The nanowire-based GaN VCSEL share much resemblance to the state of the art VCSEL described above, albeit the optical cavity looks different. The optical cavity is made of large dislocation-free hexagonal GaN nanowires, henceforth referred to as ”platelets”. A schematic view of the nanowire-based GaN VCSEL structure is presented in Figure 2.4. The platelets are put between a pair of dielectric DBRs. On the top section of the n-doped GaN platelets, a MQW structure of InGaN is grown and used as an active region. Additionally, an electron blocking layer (EBL) and a p-doped GaN layer is grown on top. The growth process of these additional layers introduces a tilted structure at the top end of the platelets. This pitch suggestively has an importance on the optical confinement and overall laser performance, and will be investigated in more detail. Figure 2.4: Schematic view of the nanowire-based GaN VCSEL. The working principle of this laser device is as follows. Electrons are injected into the device from the n-contacts and holes from the p-contacts. A current-spreading layer (ITO) is added for the purpose of leading the carriers from the p-contacts to the p-GaN. The vertical carrier confinement is enhanced due to the EBL, which ul- timately pushes the carrier recombination closer to the MQW active region. Hence, 7 2. Theory supporting carriers for stimulated emission. The optical field generated by the stim- ulated emission, however, is opposed by losses. Thus, in order to achieve lasing, the material gain must overcome these losses. As stated in the previous section, the threshold material gain (defined in Equation (2.2)) can be derived by knowing the electromagnetic field distribution within the laser cavity, thus one must solve for Maxwell’s equations. 2.2 Maxwell’s Equations and the wave equation Within the field of electromagnetism, Maxwell’s equations are used to describe the fundamental relations between some of the most important electromagnetic quanti- ties. These equations can then be used to derive the electromagnetic wave equation. Starting off with Ampère’s law ∇× ~H = ~J + ∂ ~D ∂t . (2.4) where ~H is the magnetic field intensity, ~J is the current density, and ~D is the displacement field. These quantities can be expressed in terms of the magnetic field density ~B and the electric field strength ~E, by using the constitutive relations ~B = µ0 [µr] ~H (2.5) ~D = ε0 [εr] ~E (2.6) and Ohm’s law ~J = [σ] ~E, (2.7) where the square brackets represent 3-by-3 tensors, which are used to describe gen- eral material properties. By substituting Equation (2.5), (2.6), and (2.7) into Am- père’s law Equation (2.4) and then derivating once more with respect to time, the new relation becomes ∇×  1 µ0 [µr]−1 ∂ ~B ∂t  = [σ] ∂ ~E ∂t + ε0 [εr] ∂2 ~E ∂t2 . (2.8) Using Faraday’s law, which relates the time derivative of the magnetic field density with the curl of the electric field strength ∇× ~E = −∂ ~B ∂t . (2.9) The combination of Equation (2.8) and (2.9) results in ∇× ( 1 µ0 [µr]−1∇× ~E ) + ε0 [εr] ∂2 ~E ∂t2 = − [σ] ∂ ~E ∂t , (2.10) which is the general double cross-product form of the electromagnetic wave equation. [µr]−1 is the inverse of the relative permeability tensor, [εr] is the relative permit- tivity tensor, and [σ] is the conductivity tensor. µ0 and ε0 are the scalar vacuum permeability and permittivity respectively. 8 2. Theory 2.2.1 Frequency domain representation Lasers are ideally monochromatic, hence a frequency domain representation of the wave equation is useful. Assume a harmonic time dependency of the electric field, i.e. ~E (~r, t) = ~E (~r) ejωt, (2.11) where ω is the angular frequency. Under this assumption, a time derivative operator acting on the electric field in the time domain can effectively be replaced by a multiplication in the frequency domain. Thus, ∂ ∂t ~E (~r) ejωt F→ jω ~E (~r) ejωt. (2.12) Consequently, under the above assumptions, the wave equation in Equation (2.10) reduces to ∇× [µr]−1 ( ∇× ~E ) − ω2 c2 0 ( [εr]− j [σ] ωε0 ) ~E = 0, (2.13) where c0 = √ε0µ0 is the speed of light in vacuum. 2.2.2 Physical assumptions The materials considered in this study are assumed to be isotropic. Hence, the 3-by- 3 tensors in the wave equation (2.13) reduce to scalar, spatially dependent values. The relative permittivity tensor can be expressed in terms of the material refractive index [εr] = ε = (n− jκ)2 = ( n2 − κ2 ) − 2jnκ, (2.14) where n is the refractive index, and κ is the extinction coefficient. Furthermore, the materials are assumed to be non-magnetic, thus the relative permeability tensor must be equal to unity. Also, the conductivity tensor is assumed to be equal to zero. 2.3 Optical Modes in VCSELs When applying the wave equation onto a confined structure, such as a VCSEL, the solution consists of a discrete set of field distributions called modes. Considering VCSELs, these modes are typically categorized into longitudinal and transversal modes. 2.3.1 Longitudinal Modes The optical cavity formed by the two high reflectivity distributed Bragg reflectors introduce longitudinal confinement, i.e. a confinement in the direction perpendicular to the substrate. As a consequence, the optical field becomes a superposition of forward and backward-reflected propagating waves, which form a standing wave. As a direct implication from the phase condition introduced earlier in Equation (2.1), given a separation between the top and bottom DBRs, only a discrete set of wavelengths are allowed inside the cavity: 9 2. Theory λN = 2Leffneff N , N = 1, 2, 3, ... (2.15) where Leff is the effective length of the cavity, i.e. the separation between the top and bottom mirror plus the penetration depths of the respective mirror. By derivating the above relation with respect to N , and then let |∆N | = 1, it can be shown that the separation between the longitudinal modes is given by |∆λ| = λ2 2neffL [ 1− λ neff ( ∂neff ∂λ )] . (2.16) Since the cavity length in a VCSEL is typically very short, the separation between the longitudinal modes is large. Hence, in combination with the band-limiting effects of the DBRs, the VCSEL usually operates with only one single longitudinal mode [9]. 2.3.2 Transversal Modes The transversal optical confinement (radial direction) in a VCSEL requires some radial variation in local refractive index, such that the local refractive index is high for small r and low for large r. If these conditions are met, the optical cavity will exhibit wave-guiding properties. In the case of the platelet-based VCSEL, a radial index step is a natural consequence of the small finite width of the platelet. In order to adjust the transversal confinement, the platelet is surrounded by a cladding layer. By further assuming that the cavity is cylindrical and has infinite length (or equivalently, open ends), the platelet structure resembles the waveguide of a step- index optical fiber1[9]. The field distributions of such step-index optical fiber are well understood, and a derivation of the corresponding modes is presented in Appendix A. The six lowest order optical fiber modes are presented in Figure 2.5. The field distributions presented in the figure have been generated in the simulation software COMSOL Multiphysics®. There are three types of modes. Some modes are fully transversal with no longitudinal component. Those are called transversal electric (TE0p) if Ez = 0 or transversal magnetic (TM0p) if Hz = 0. Modes of the third type have non-zero components even in the longitudinal direction, thus they are called hybrid modes and are abbreviated as HEmp. The index p is the solution index and m is the azimuthal wavenumber. It can be shown that the first order mode is HE11, which has an azimuthal wavenumber m = 1. This knowledge will be important later on when setting up the platelet simulations in COMSOL Multiphysics®. 1The same argument goes for other conventional VCSEL structures as well. 10 2. Theory HE11 0 0.2 0.4 0.6 0.8 1 HE11 0 0.2 0.4 0.6 0.8 1 TE01 0 0.2 0.4 0.6 0.8 1 HE21 0 0.2 0.4 0.6 0.8 1 HE21 0 0.2 0.4 0.6 0.8 1 TM01 0 0.2 0.4 0.6 0.8 1 Lowest order mode First set of higher order modes Figure 2.5: First six transversal modes of a cylindric step-index waveguide with radius Rpl = 450 nm in an infinitely long GaN-platelet surrounded by air. The refractive indecies are ncore = 2.48 and nclad = 1.5, and the cladding thickness is 100 nm. The colormap indicates the normalized norm of the electric field, and the white arrows indicate direction and strength of in-plane electrical field components. 2.3.3 Cylindrical geometry approximation The GaN platelets have a hexagonal cross section, which complicates the modelling and a 3D-model is required to fully represent the structure. Thus, in order to save computational time and effort, the hexagonal platelets have been approximated as cylinders. The rotational symmetry gained from this approximation ultimately reduces the aforementioned 3D-problem to an axial-symmetric 2D problem. The effects of changing geometry from hexagonal to circular cross sections have been investigated by Henneghien et al. (2009)[10]. They have shown by the means of finite difference time domain (FDTD) simulations that nanowires with hexagonal cross sections share many optical properties to that of circular structures. They also concluded that the main differences derive from a difference in cross-sectional area. If the cross sectional-area is kept the same, however, the similarities are very strong. Therefore, all simulations performed in this thesis have been run under the cylindrical approximation. These results are valid for cylindrical cavities, but the results from Hennighien et al. show implications that the field distributions and trends derived in this thesis can be translated to hexagonal structures. To translate the cylindrical geometry to a comparable hexagonal geometry the following equation can be used: 3 √ 3 2 a2 = πR2 pl ⇒ a ≈ 1.1Rpl (2.17) where a is the side length of the equilateral hexagon, and Rpl is the radius of the 11 2. Theory HE11 0 0.2 0.4 0.6 0.8 1 HE11 0 0.2 0.4 0.6 0.8 1 TE01 0 0.2 0.4 0.6 0.8 1 HE21 0 0.2 0.4 0.6 0.8 1 HE21 0 0.2 0.4 0.6 0.8 1 TM01 0 0.2 0.4 0.6 0.8 1 Lowest order mode First set of higher order modes Figure 2.6: First six transversal modes of a hexagonal step-index waveguide with side length a ≈ 1.1Rpl in an infinitely long GaN-platelet surrounded by air. Rpl is set to 450 nm, and the cladding thickness is 100 nm. The refractive indecies are ncore = 2.48 and nclad = 1.5. The colormap indicates the normalized norm of the electric field, and the white arrows indicate direction and strength of in-plane electrical field components. circular platelet. Thus in order for the two geometries to have the same cross- sectional area a ≈ 1.1Rpl. The first six transversal modes of a cylindrical cavity of radius Rpl = 450 nm and cladding (100 nm) in Figure 2.5 can be compared to the first six transversal modes of a hexagonal cross sectional cavity with side length a ≈ 1.1Rpl in Figure 2.6. The refractive indices of the core and cladding are in both cases set to ncore = 2.48 and nclad = 1.5 respectively. Both structures have the same cross-sectional area, and the mode shapes are very similar, especially the fundamental mode HE11. The fundamental mode is the main focus of this thesis, since it usually experiences the lowest losses. Since the mode distributions are very similar for this mode for cylindrical and hexagonal waveguides, it is a clear indication that the results and trends explored in the cylindrical case, can be translated into the hexagonal case. 2.4 The Finite Element Frequency Domain Method The finite element method (FEM) is a numerical method designed to solve partial differential equations. One of its key features is that the method can be used to model most arbitrary and complex shapes and structures while still providing accu- rate results. Furthermore, the matrices associated with the method are very sparse, which leads to reduced memory consumption. Hence, the FEM is well established 12 2. Theory and used in many applications[11]. The finite element frequency domain (FEFD) method is developed to solve the time harmonic electromagnetic wave equation by numerical means. The basic idea of the FEM and FEFD is to discretize an infinite dimensional problem into a finite dimensional problem. This is done in the spirit of Galerkin by the introduction of the so called weak form of the differential equation. First, it is assumed that the electrical field can be described as a linear combination of test functions ~Nµ ~E = ∑ µ Eµ ~Nµ. (2.18) The test functions form a basis set. The wave equation is then multiplied with the test function ~Nν and integrated over the entire space. This is called the Galerkin approach, and it will after some partial integration result in the following so called weak form.∫ V 1 µr ( ∇× ~E ) ( ∇× ~Nν ) dV − ∫ V ω2 c2 0 εr ~E · ~Nν dV = − ∮ ∂V 1 µr n̂× ( ∇× ~E ) · ~Nν ds. (2.19) where, V is the volume of interest and ∂V is its boundary. Here, we have made the physical assumptions described in section 2.2.2, i.e. isotropicity and no conductivity. By expanding ~E in terms of the test functions ~Nµ described in Equation (2.18), the weak form is reduced to a matrix eigenvalue problem MµνEµ + λ̃2KµνEµ = fν , ν = 1, 2, 3, ..., Ndof (2.20) where Mµν = ∫ V 1 µr ( ∇× ~Nµ ) ( ∇× ~Nν ) dV (2.21) Kµν = ∫ V εr c2 0 ~Nµ · ~Nν dV (2.22) fν = − ∮ ∂V 1 µr n̂× ( ∇× ~E ) · ~Nν ds (2.23) Due to the equation’s form, Mµν is the called the mass matrix, Kµν is the stiffness matrix and fν is the load vector. This is a set of Ndof equations, where each ν correspond to one degree of freedom. By selecting basis functions with local support, i.e. they are non-zero only in a small geometrical region, the overlap between basis functions become small. Since, the integrals become non-zero only when there is an overlap between two basis functions, the system matrices become sparse, which is one of the key benefits of the FEM. The eigenvalue λ̃ is the angular frequency of the system. By allowing complex eigenvalues, λ̃ = jω + δ, (2.24) the imaginary part will correspond to oscillations, and the real part will represent attenuation (if δ < 0) since ~E (~r, t) = ~E (~r) eλ̃t = ~E (~r) e(jω+δ)t = ~E (~r) ejωt︸︷︷︸ Oscillation · eδt︸︷︷︸ Attenuation . (2.25) 13 2. Theory Hence, by solving this eigenvalue problem, the eigenvalues will provide information about the resonance frequencies and also the attenuation (damping). Provided with both the damping and the field distributions, FEFD can be used to determine the threshold material gain of the platelet-based VCSEL. In order to obtain accurate results, the structure must be properly meshed. But choosing size and shape of the individual elements is an art of its own. Hence, the simulation of the laser structure has been performed using the FEFD technique in the software program COMSOL Multiphysics®, which comes with meshing algorithms designed to optimize the element features. 14 3 Methods This section describes the programs and custom algorithms used and developed for the purpose of examining the optical properties of the platelet-based VCSEL. The first part introduces COMSOL Multiphysics® and how to set up the program to run optical simulations. The second part describes the use of LiveLinkTM for MAT- LAB®, which makes it possible to control COMSOL® simulations with MATLAB® commands. The third part presents the layout and principles of the main program and the last part describes the post-processing algorithms used to extract the results. 3.1 COMSOL Multiphysics® Simulation of the laser structure has been performed using finite-element frequency domain technique in COMSOL Multiphysics®. A simulation in COMSOL® involves six steps 1. Define geometries 2. Implement physics 3. Generate a mesh 4. Configure eigenvalue solver 5. Run simulation 6. Post-process simulation output The following sections will go through and explain these steps in more detail. 3.1.1 Geometries of interest As discussed earlier, the hexagonal platelet-based VCSEL-structure can be treated as a rotationally symmetric cylinder. Hence, the 3D-problem is reduced to a 2D- axisymmetric problem. We here on out work in cylindrical coordinates centered at the origin of the platelet. A geometric sketch of the platelet-based VCSEL is presented in Figure 3.1. 15 3. Methods Figure 3.1: Geometry of platelet-based GaN VCSEL. Once the dimension of the problem is set, the geometry of the VCSEL-structure must be defined. In COMSOL®, the geometry can either be imported from an ex- ternal computer aided design (CAD) tool or be constructed from scratch directly in the COMSOL Desktop® window. In this study, the geometry was constructed from scratch. The radius and thickness of each region is defined from a set of global parameters. These are all accessible and can be adjusted when needed. The param- eters are categorized in three groups: geometrical, radial, and vertical parameters. A list of these parameters together with their default values are presented in Tables 3.1, 3.2, and in 3.3. Table 3.1: List of Geometrical Parameters Parameter Default value Description θ 58◦ Tilt angle on platelet top section nbrOfTopDBR 12 Number of mirror pairs in top DBR nbrOfBottomDBR 10 Number of mirror pairs in bottom DBR isFlatDBR true Flat or concave top DBR δactive 0 nm MQW position displacement from default 16 3. Methods Table 3.2: List of Vertical Parameters Parameter Default value [nm] Description dsapphire 420 Sapphire substrate thickness dDBR,low 71.4 Low index DBR layer thickness dDBR,high 49.8 High index DBR layer thickness dspacing 31 Spacing layer thickness dITO 30 ITO (indium tin-oxide) thickness dpGaN2 140 p-GaN length (above EBL) dEBL 20 EBL thickness dpGaN1 13 p-GaN length (below EBL) dQWw 5 Quantum well well thickness dQWb 6 Quantum well barrier thickness dnGaNp 27 n-GaN length (pitched only) dnGaNv 770 n-GaN length (vertical only) dnGaN dnGaNv + dnGaNp n-GaN length dmask 30 Si3N4-mask thickness dsubstrate 100 GaN substrate thickness dair 420 Air region thickness dpl 1009 Platelet height Table 3.3: List of Radial Parameters Parameter Default value [nm] Description Rsimbox 2000 Simulation Box Size RPML λ0/2 Thickness of PML Rpl 450 Platelet radius Rclad 100 Cladding layer thickness Rmask 25 Mask hole radius RplTop Rpl − (dpl − dnGaNv) ∗ cot ( θ π 180 ) Platelet radius at top end Rembedding Rsimbox −Rpl −Rclad Embedding layer width 3.1.2 Implementing Physics The physics package of interest is ”Electromagnetic Waves, Frequency Domain (ewfd)”, and it has been used to solve the wave equation for time-harmonic fields. Starting from the Equation (2.13) and taken into account the assumptions stated in Section 2.2.2, the wave equation reduces to ∇× ( ∇× ~E (~r) ) − ω2 (n− jκ)2 c2 0 ~E (~r) = 0 (3.1) where ~E is the electric field, ω is the angular frequency, n is the refractive index, κ is the extinction coefficient, and c0 is the vacuum speed of light. The material properties can be set to be regionally dependent, which allows for optical simulations in complex structures such as the platelet-based VCSEL. The 17 3. Methods materials in the simulation are represented by specifying a refractive index and an extinction coefficient. A list of the material parameters used in the simulations are presented in Table 3.4. Table 3.4: List of Material Parameters Parameter Default value Description nQWb 2.68 Refractive index of In0.10Ga0.90N in QW barrier nQWw 3.00 Refractive index of In0.03Ga0.97N in QW well nEBL 2.42 Refractive index of EBL nGaN 2.49 Refractive index of GaN nITO 2.10 Refractive index of ITO nair 1.00 Refractive index of air nembedding 1.00 Refractive index of embedding material nmask 2.06 Refractive index of mask nsapphire 1.78 Refractive index of sapphire substrate nClad 1.50 Refractive index of cladding nSiO2 1.47 Refractive index of SiO2 in DBR1 nHfO2 2.11 Refractive index of HfO2 in DBR1 The material absorption is related to the imaginary part of the refractive index, i.e. the extinction coefficient, i.e. κ = α λ0 4π (3.2) where α is the absorption coefficient and λ0 is the vacuum wavelength. Consequently, absorption in the materials can be added to the simulation through defining κ as per above. Also, absorption can be switched on or off by simply multiplying the extinction coefficient by a boolean parameter isAbsOn. Per default the absorption is turned off. 3.1.2.1 Boundary Conditions Consider a system with open boundaries. That could be for instance a system with some type of source irradiating into infinite space, or a planar laser cavity with planes that seemingly extend into infinity. When performing electromagnetic simulations on such systems with open boundaries, it is necessary to truncate the simulation region, due to practical reasons. This truncation, however, introduces unwanted numerical reflections at the boundaries which interfere with the field distribution causing unreliable results [12]. This implies that in order to suppress the artifacts caused by truncation, the boundaries must be both absorbant and non-reflecting. Simple absorbing boundary conditions, such as Scattering Boundary Condition (SBC) in COMSOL®, perform well and are practically non-reflecting and fully ab- sorbent for waves with normal incidence. Albeit, for waves with oblique incidence, 1Photonics Laboratory, Department of Microtechnology and Nanoscience, Chalmers University of Technology 18 3. Methods the SBC breaks down and artificial reflections occur [13]. The SBC works well in 1D and in cases when the waves propagate normal to the boundaries. For the nanowire- based GaN VCSEL structure, early simulations performed within the scope of this thesis, imply that the cavity suffers from scattering losses, thus the waves will most likely hit the boundaries with oblique incidence. Therefore, another approach for dealing with artificial reflections must be used. 3.1.2.2 Perfectly Matched Layers (PML) A more promising approach is that of using Perfectly Matched Layers (PML). A PML is not a ”boundary” condition but rather a layer with absorbing and non- reflective properties, see Figure 3.2. The PML used in COMSOL® is based on Coordinate Stretching, which is defined as follows. First, let’s perform an analytical continuation from real to complex coordinates inside the PML region x→ { x, x < x0 x+ F (x), x0 < x < xedge (3.3) where x0 marks the interface between the physical domain and the PML, and xedge is the edge of the simulation box. In COMSOL® F (x) is defined as[14] F (x) = λtf(ξ)−∆wξ, ξ ∆= x− x0 ∆w (3.4) where λt is the ”typical wavelength”, and ∆w = xedge − x0 is the thickness of the PML. For eigenvalue studies the typical wavelength is set to equal the PML thickness by default, i.e. λt = ∆w. f(ξ) = f ′(ξ) − jf ′′(ξ) is the complex PML stretching function expressed in a normalized coordinate ξ, which is equal to 0 at x = x0 and 1 at the edge of the simulation box x = xedge. x x0 xedge ∆w Physical Domain PML Figure 3.2: Schematic illustration of a 1D PML Consider a propagating planewave with wavevector kx inside the PML. By im- plementing the coordinate stretching Equation (3.3) and (3.4), the plane wave can be factorized into three parts that includes the incident plane wave at the PML border, one oscillatory factor and one attenuation factor. e−jkxx → e−jkxxe−jkxλtf(ξ)e+jkx∆wξ = e−jkxxe−jkxλtf(ξ)e+jkx(x− x0) = e−jkxx0e−jkxλtf(ξ) = e−jkxx0︸ ︷︷ ︸ Incident plane wave e−jkxλtf ′(ξ)︸ ︷︷ ︸ Oscillations e−kxλtf ′′(ξ)︸ ︷︷ ︸ Attenuation (3.5) 19 3. Methods It is clear that the attenuation is scaled by the imaginary part of the stretching function, and the real part creates additional oscillations inside the PML. Hence in order to damp the wave, a large imaginary part of the stretching function is prefer- able. However, if the wavevector has an imaginary part, i.e. evanescent waves are present, the issue gets more complicated. For evanescent waves, a large imaginary part of the stretching function cause faster oscillations and instead a large real part of the stretching function is preferable. Consequently, one has to find a balance with not too high nor too low real and imaginary parts. The default PML stretching function in COMSOL® is of polynomial type f(ξ) = f ′(ξ)− jf ′′(ξ) = sξΠ(1− j) (3.6) where s is called the PML scaling factor, and Π is the curvature factor. Per default, s = 1 and Π = 1. A high s results in faster decaying fields, but it also cause more oscillating solutions. In conclusion, the scaling factor should be kept close to s = 1. Increased curvature of the PML stretching function requires a denser grid. Otherwise the solution can become subject of under-sampling, which may cause divergences at the PML boundary. The curvature parameter has therefore been set to Π = 1. Consequently, the default values have turned out to work well. The thickness of the PML layer ∆ω = RPML is chosen to be half wavelength as a rule of thumb stated in [12]. By doing so the PML serves as an absorbing and non-reflective boundary. It is worth noting that the implementation of PML boundaries introduces ar- tificial field solutions amongst the real modes. Their appearance is characterized by having strong field components close to the simulation box boundary, and also having strong oscillations. These solutions should and have been ignored. 3.1.3 Meshing A good quality mesh is crucial for obtaining reliable solutions. COMSOL Multi- physics® has built-in meshing algorithms that generates optimized meshes depend- ing on the geometry. It also provides the option to generate a physics-defined meshes that adapt to the local material properties. But in order to be able to control the mesh properties, the mesh used in this work has been generated with customized directives. In order to properly resolve the propagating waves in an optical simulation, the maximum element size, h, must be less than a fraction of the local wavelength, that is h ∼ 1 M λ0 nlocal(~r) , M = 8, 10, 12 (3.7) where nlocal(~r) is the refractive index of the local material. This is an empirical estimate for a sufficient element size to yield accurate results. Too small elements, however, might result in very long computation time though. In 2D, the elements can have many different shapes: triangular, quadrilateral, polygonal. Generally, a triangular mesh is preferable for many applications. But for structures which exhibits lots of symmetries, for instance 2D-axisymmetry, one should consider using quadrilateral meshing along the symmetry axis [15]. 20 3. Methods Early meshing attempts with entirely triangular mesh occasionally resulted in nonphysical solutions. Those had a singularity, i.e. a sudden localized high inten- sity peak situated close to the symmetry axis. A smaller grid size did not solve the problem. But by following the advice in[15], and made sure that the edge elements along the symmetry axis were of quadrilateral shape, the singularities were sup- pressed. Therefore, the mesh in the subsequent simulations were of a hybrid type with triangular elements in the curved structures and quadrilateral in the rest. 3.1.4 Solver Configuration The eigenvalue problem associated with the finite element method has a large num- ber of solutions, most of which are irrelevant for specific studies. Thus, there are two search methods in COMSOL® that help narrow down the search for relevant solutions. There is a "regional" search method, and a "manual" search method[14]. The regional search finds all eigenvalues within a user-specified region in the com- plex plane. This method has an option to perform a consistency check, which makes sure that the solver doesn’t unintentionally miss any eigenvalues inside the region in the initial run. This is a neat feature, however, the consistency check requires additional computational time. Furthermore, the user-specified region might include more eigenvalues than what is considered relevant for the study, thus, it will most likely also result in an excessive use of computational power and time. Ultimately, the manual search option was preferred above the regional search method. In the manual search method, the user specifies a reference point and a fix number of requested eigenvalues, Neigs. The solver then provides only the Neigs num- ber of eigenvalues that are the closest to the reference point. This method requires no consistency check, and the number of found eigenvalues are fixed. Consequently, the computational time is kept low. Both search methods above require an initial reference point. It is therefore a necessity to at least somewhat know in which region of the wavelength spectrum one expects to find relevant solutions. Since the structure and its materials are known, it is possible to predict an approximate wavelength by approximating the structure as a simple Fabry-Perot cavity, which consists of one uniform material put between a pair of perfectly reflecting mirrors. By implementing the phase condition stated in Equation (2.1), the resonance wavelength can be estimated from a few known parameters. The distance between the to DBRs is L ∼ 1200 nm. That includes the platelet, the GaN substrate, the Si3N4-mask, the ITO, and the spacing layer. Furthermore, from the refractive indices of the core materials, it is reasonable to assume that the effective refractive index should lay between 2.09 (HfO2) and 3.00 (In0.03Ga0.97N). Most of the field propagates within GaN, thus the effective refractive index should lay very close to 2.5. In addition, since the low index materials in the cavity are more abundant than the high index materials, the effective index is slightly less, and it is therefore assumed to be around 2.45. A GaN VCSEL structure with strong similarities to the platelet-based GaN VCSEL, has been studied by Å. Haglund et. al. in 2016 [1]. The cavity shown in the paper is based on GaN with InGaN MQW as active region, and is put between a pair 21 3. Methods of dielectric DBRs. With a cavity length of about 1200 nm the fields presented in the paper exhibited 14 optical field anti-nodes between the two DBRs. Consequently, by also assuming 14 optical field anti-nodes for the platelet-based GaN VCSEL, the resonance wavelength is expected to be approximately λref ≈ 2neffL N = 2× 2.45× 1200 nm 14 ≈ 420 nm. (3.8) With the above prediction the reference eigenfrequency throughout the proceed- ing simulations has been set to fref = c0 λref = 299′792′458 m/s 420 nm ≈ 7.1379× 1014 Hz. (3.9) The resonant wavelength will slowly drift as the cavity structure is varied. Hence, to make sure that the resonant wavelength is found among the fixed number of solutions even for a significant change of cavity structure, Neigs is set somewhat higher than default. Typically, Neigs = 20. 3.2 LivelinkTM for MATLAB® The COMSOL Desktop® software has a graphical user interface that works well. As discussed above, COMSOL® supports CAD tool for generating geometries, physics- packages designed for optical simulations, and a in-built meshing algorithm that generates high quality meshes. COMSOL® also has post-processing tools for visu- alizing the simulation results. But when it comes to more complex post-processing that involves combinations of data from different simulation runs, COMSOL Desk- top® is a bit convoluted. For example if results are exported to an external software for data processing, such as MATLAB® or Excel® one has to manually export the data after each simulation step. Through looking for other alternatives, LivelinkTM for MATLAB® turned out to be a good alternative. With LivelinkTM for MATLAB® it is possible to connect and control COMSOL Multiphysics® remotely in a MATLAB® scripting environ- ment, which allows more complex set-ups and post-processing. Hence, LivelinkTM for MATLAB® lays the foundation for the program that calculates the threshold material gain. 3.3 The Main Program The main program has three important purposes. It must be able to establish a connection to and communicate with the COMSOL Server®. It has to be able to perform parameter sweeps for various structural and material parameters. And it has to be able to extract and save relevant data from the simulation suitable for post-processing. A flow chart that describes the main program is presented in Figure 3.3. The script starts by establishing a connection to the COMSOL Server®. Next, geomet- rical, material and solution configurative parameters are specified by the user. The 22 3. Methods user also setups sweep parameters together with a specified sweep interval. In the current version, up to two sweep parameters can be selected simultaneously, which allows for 2D-mapping and 3D-visualization. Outer Loop Parameter Sweep Gain Sweep Connect to COMSOL Server® Setup simulation • Filenames • Custom parameters • Sweep intervals Set outer loop parameter(s) Set active gain gactive Run simulation Done? Extract data Eigenvalues λ̃ • Field distributions • Mesh • VCS profiles Done? Finish Yes Yes No No Figure 3.3: Flow chart of main program main.m. Once the simulation has been set up, the script enters a nestled for-loop that sweeps over the sweep parameters. One or two outer loops to sweep over sweep pa- rameters, and one inner loop that extracts data necessary to determine the threshold material gain. This inner loop is called a gain sweep and its purpose and function is described in more detail in the next section. Once the gain sweep is finished rele- vant data is extracted from the COMSOL® model. That includes eigenvalues, field distributions, meshes and vertical cross-section (VCS) profiles. The eigenvalues are saved in .dat file format, and field distributions for all found solutions together with the associated mesh are saved in .fig file format to allow for later analysis. The field data along the z-axis is saved separately in a MATLAB® struct that can be used to form VCS profiles. The MATLAB® struct also contains all descriptive information 23 3. Methods about the simulation configuration including all parameters, the sweep intervals and the simulation time. 3.3.1 Gain Sweep Algorithm After importing the VCSEL structure into COMSOL®, the wave equation is numer- ically solved as an eigenvalue problem. The output eigenvalue is a complex angular eigenfrequency λ̃ = −jω + δ (3.10) where ω is the angular eigenfrequency, and δ is the damping of the solution. The damping is directly related to losses in the device. Hence, the damping must also be related to the threshold material gain. The damping is expressed in [rad/s] while the threshold material gain has units of [1/cm]. Therefore, it is necessary to find a way to translate between the two quantities. A straight forward way to estimate the threshold material gain is to apply a gain to the active region, and then successively increase it until the losses, i.e. the damp- ing term, vanishes. This is potentially a time consuming method, since the entire eigenvalue problem must be solved multiple times. However, the relation between the applied gain and the damping term turned out to be simple and predictable. Consequently, this method became the main approach to determine the threshold material gain. An initial gain sweep was performed for a structure with platelet radius Rpl = 450 nm, cladding thickness of Rclad = 100 nm, and cladding refractive index of nclad = 1.5. The platelet was here assumed to be fully cylindrical with no pitch in the top section, and the absorptive properties of the materials were turned off. The simulation box-radius was set to RSimbox = 3000 nm, and the initial wavelength guess was set to λref = 411, 76 nm, which had been the output wavelength from a previous run. Lastly, the perfectly matched layer had a thickness of λref/2, and the scaling and curvature were set to unity. Figure 3.4 shows the damping term as it decreases with a higher applied gain. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 10 4 -5 0 5 10 15 20 10 11 Figure 3.4: Damping term as a function of applied gain to the active region, including a least square fit to the simulated δ-points. 24 3. Methods The relation between the applied gain and the damping term is linear and can be expressed as δ = kg + δ0. (3.11) Thus, the threshold gain can be determined from this linear relation by setting δ = 0 and solve for the applied gain gth = −δ0 k (3.12) Consequently, in order to determine the threshold gain, all that is needed is an estimation of the y-intercept δ0 and the slope k. This estimation was done using a least square linear regression algorithm, which is based on minimizing the sum of the square of the residual Q. By assuming the linear relation (3.11) between the response δi and the predictor gi, Q is defined as Q = ∑ i (δi − (δ0 + kgi))2 . (3.13) By minimizing Q with respect to δ0 and k, these coefficients can be determined as follows k = ∑ i (gi − ḡ) ( δi − δ̄ ) ∑ i (gi − ḡ)2 (3.14) δ0 = δ̄ − kḡ (3.15) where δ̄ and ḡ represent the mean value of the sample points (gi, δi), i.e. δ̄ ∆= 1 Ngs Ngs∑ i=1 δi (3.16) ḡ ∆= 1 Ngs Ngs∑ i=1 gi. (3.17) In order to obtain the threshold material gain for one parameter configuration, there is need for a gain sweep. First a gain g1 is applied to the active region via the extinction coefficient κ = −gλref4π . (3.18) The simulation then runs and solves the eigenvalue problem and outputs Neigs num- ber of eigenvalues. These are then stored. The applied gain is changed to g2 > g1 and the simulation runs again. This procedure is repeated for a few steps. The eigenvalues from the different gain sweep steps are then sent to the post-processing algorithm, which sorts and extracts the applied gain vs damping relation. Finally, the threshold material gain can be determined from those curves by the Least Square Linear Regression Algorithm described above. Due to the linear relation only a few sample points is required for achieving consistent results. Hence, there is no need to step-wise increase the applied gain until the damping term vanishes, thus the simulation time is manageable. 25 3. Methods 3.3.1.1 Threshold correction Equation (3.18) describes the gain’s wavelength dependence. Both gain and damp- ing depends on the absolute wavelength. Therefore, the deviation between the wavelength guess λref and the actual wavelength λ0, must be accounted for. The difference in apparent and actual gain due to wavelength mismatch is g̃ = −4π λ0 κ = g ( λref λ0 ) . (3.19) The damping factor δ is directly derived from the eigenvalue solver, so there is no mismatch in δ. However, the mismatch in apparent and actual gain will affect the coefficients in the gain-damping relation δ̃ = δ = kg + δ0 = k̃g̃ + δ̃0. (3.20) Since g = 0 implies g̃ = 0, the y-intercept is consistent in both apparent and actual cases, i.e. δ̃0 = δ0. The slope, however, scales the same way as the applied gain. Consequently, the actual threshold material gain becomes g̃th = −δ0 k̃ = − ( λref λ0 ) δ0 k = gth ( λref λ0 ) (3.21) In conclusion, in order to obtain the actual threshold gain one must multiply the apparent threshold with the fraction between the wavelength guess and the actual wavelength. 3.4 Post-Processing Algorithms The post-processing script has two main objectives: to extract the threshold ma- terial gain using least square linear regression, and to plot field distributions. The threshold material gain can be extracted directly from the eigenvalue output of the main program. Unfortunately, the output can sometimes become slightly unsorted, which requires some tedious sorting procedures to fix. At each gain sweep step Neigs eigenvalues are output to the data file. The eigen- values are put in rows and the eigenfrequency solver in COMSOL® automatically sorts the eigenvalues in ascending order. For each step, a new column of eigenvalue output is printed. If assumed that the wavelength for each mode is fixed upon gain stepping, this spreadsheet format makes it straight forward to determine the thresh- old material gain by simply read the real part of the eigenvalue from each gain sweep step (column) and perform the extraction method described in Section 3.3.1. In addition to the wavelength drifting deriving from geometrical parameter changes, it has been noticed that the wavelength also drifts slightly when different gain is applied to the active region. This phenomena is the cause of an unwanted ef- fect that henceforth is referred to as "list-shifting". There are two types of list-shifts. First, upon gain sweep stepping, one eigenvalue inside but close to the perimeter of the solver search region might end up being replaced by a new eigenvalue that previously was located outside the search region, see Figure 3.5. This cause the 26 3. Methods =⇒ fref fref • • • • • • • •◦ ◦ ◦ ◦ g1 g2 Figure 3.5: The origin of type I list shifts. Eigenvalues drift due to a change in applied gain. This can eventually cause two eigenvalues close to the perimeter of the search region to switch places, which causes a shift in the list of eigenvalues, hence the name. entire row of eigenvalues from one gain sweep step to shift in relation to the other steps. Second, it may occur that two separate eigenvalues with similar wavelength switch places upon gain sweep stepping, which cause a local shift. In order to handle the list-shifts, a group and sort algorithm has been imple- mented. But to avoid unnecessary sorting on data sets not contaminated by list- shifts, a list-shift detection algorithm has been implemented. 3.4.1 List-shift detection First, the damping is extracted row by row which forms a set of damping-gain curves. If no list-shift has occurred, those are linear curves and the extraction of threshold material gain is straight forward. However, if a list-shift has occurred, the lines are entangled, see Figure 3.6. The figure also imply that a clear characteristic of a list-shift is a sudden change in slope, thus the local point-wise slope has been used as tool for detection of list- shifts. The local point-wise slope has been defined using a forward finite difference scheme kli = δli+1 − δli gi+1 − gi , i = 1, 2, ..., NGS − 1, l = 1, 2, ..., Neigs, (3.22) where NGS is the number of gain sweep steps, i.e. columns, and Neigs is the number of requested eigenvalues per gain sweep step, i.e. rows. Furthermore, a row-by-row standard deviation of the point-wise slopes is calculated. Sudden changes in point- wise slope due to a list-shift result in a large standard deviation, several orders of magnitude larger than if no list-shift had occurred. Hence, by calculating the maximum of the standard deviations for each outer loop parameter step and then compare it to a tolerance level, it is possible to detect the potential presence of list-shifts, i.e. list-shifts are present if max l { σl ( kli )} > tolerance. (3.23) 27 3. Methods 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 10 5 0 2 4 6 8 10 12 14 10 13 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 10 5 0 2 4 6 8 10 12 14 10 13 Figure 3.6: Illustration of typical damping-gain curves. Top: with list-shifts. Bottom: no list-shifts. 3.4.2 List-shift processing After the data sets which have been affected by list-shifts have been detected, the unsorted data will be processed. The eigenvalues are first rearranged from the de- fault 2D spreadsheet format into a 1D array sorted with respect to ascending angular frequency. From the assumption that the eigenvalue shift due to gain sweeping is sig- nificantly smaller than the eigenvalue separation between different mode solutions, the eigenvalues are grouped by angular frequency, i.e. eigenvalues with similar angu- lar frequencies are grouped together. Some groups will have less than NGS members, which is an indication of type I list-shifts. The unfilled groups will be excluded from the set. The remaining groups might suffer from type II list-shifts, thus a simple sorting and grouping with respect to only angular frequency won’t be enough. Consequently, a second sorting and grouping is needed where also the local slope is evaluated. As pointed out above, the local slope k from eigenvalues corresponding to the same mode solution should be constant. Thus if the evaluated local slopes for two sub- sequent eigenvalues differs more than a set tolerance, the eigenvalues correspond to different modes and are therefore separated into different groups. After this second grouping procedure, all remaining unfilled groups and unpaired eigenvalues are sent 28 3. Methods and treated in a third grouping. In the third grouping, the remaining eigenvalues are first again sorted with respect to ascending angular frequency. Next, the first eigen- value is taken as starting point and is paired with the eigenvalues that are closest in absolute value. Eigenvalues that have been paired are removed from the list. If the group is filled, a new starting point is selected from the remaining list, and so on. Since all groups which suffered from type I list-shifting were excluded after the first grouping, all groups from the second and third grouping will be filled. Consequently, in the end, all remaining groups will contain sorted eigenvalues with corresponding applied gain from which the threshold material gain can be determined. 29 3. Methods 30 4 Results This section presents the results from the parameter sweeps performed in this the- sis. The effect on the threshold material gain due to the platelet radius, cladding thickness and cladding refractive index, has been studied, as well as for two types top DBRs and for small displacements of the active region. In the last section, an error estimation study is presented in order to validate the results. 4.1 Platelet Radius From the analytical derivation of the transversal modes of a step-index optical fiber in Appendix A, it is implied that the platelet radius is a very important parameter that directly affects the wave-guiding properties of the laser cavity. In order to evaluate the effect of the platelet radius on the threshold material gain of the platelet- based GaN VCSEL, a platelet radius parameter sweep has been performed. The platelet radius was swept in the interval 300 nm ≤ Rpl ≤ 1500 nm. Furthermore, 400 600 800 1000 1200 1400 0 0.5 1 1.5 2 2.5 3 10 5 900 1000 1100 1200 1300 1400 1500 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 Figure 4.1: Threshold material gain as a function of platelet radius for a pitched top and a straight top. 31 4. Results the influence on the threshold material gain due to the pitched top of the platelet has also been studied and compared to a straight top case. The realistic pitched top has an tilt angle θ = 58 ◦, and the straight top corresponds to θ = 90 ◦. For the whole sweep the number of DBR mirror pairs were 15 for the top DBR and 10 for the bottom DBR. The simulation box radius was set to Rsimbox = 2000 nm, and the PML thickness was set to half wavelength RPML = λref/2 = 210 nm. The mesh was a triangular-quadrilateral hybrid with maximum element size h ∼ 1/12× λref/nlocal. The results from the platelet radius sweep can be seen in Figure 4.1, and the corresponding field distributions for some interesting points are shown in Figures 4.2 through 4.4. Figure 4.2: Electric field distributions for platelet radius Rpl = 450 nm. Pitched top (left) and straight top (right). 32 4. Results Figure 4.3: Electric field distributions for platelet radius Rpl = 800 nm. Pitched top (left) and straight top (right). The radius of the platelet has shown significant effect on the threshold material gain. Also the pitched top has consistently higher threshold compared to the results from the straight top. A possible explanation for this behaviour can be found by looking at the field distributions. For small radii, the optical cavity has dimensions that are comparable to the wavelength of the light, thus diffraction effects are significant, as can be seen by looking at the field distributions in Figure 4.2. The confined cavity causes wide scattering angles, hence the losses are high and the thresholds are far beyond any practical levels. Furthermore, the pitched top is shown to suffer more from diffrac- tion losses compared to the straight top case, which is concluded from the clearly visible rays inside the DBRs, hence also the difference in threshold material gain. 33 4. Results Figure 4.4: Electric field distributions for platelet radius Rpl = 1300 nm. Pitched top (left) and straight top (right). At a radius of about Rpl ∼ 800 nm the dimensions of the platelet are about four times larger than the wavelength and both diffraction and scattering are reduced, see Figure 4.3. The field start to resemble standing waves. The thresholds are still very high, but they are descending as the platelet radius is further increased. At about Rpl ∼ 1300 nm, the threshold has reached levels that could be provided by the QWs. By further increasing the platelet radius, the threshold is reduced, which agrees with our expectations. The trend can be explained as the cavity radius tends to infinity, the NW laser more and more resembles an ordinary laser from planar materials. In addition, the pitched top becomes less and less significant as the platelet radius is increased, which is also indicated by the similarities in the field distributions in Figure 4.4. 34 4. Results The platelet radius simulations suggests that large radius platelets are generally preferable over platelets with small radius, due to generally lower threshold material gains. These extra-wide platelets are however more challenging to grow, thus it is of interest to try finding complementary design alternatives that reduce the threshold material gain and that do not require too large platelet radii. 4.2 Cladding thickness and refractive index A simulation study of the cladding layer thickness and refractive index has been performed. The platelet radius was pinned at Rpl ∼ 1300 nm, and only realistically tilted structures were considered. The cladding thickness, RClad, was swept between 100 nm and 300 nm. The choice of cladding material determines the refractive index, which here was varied in between 1.5 and 2.6. SiO2 provide a refractive index at about 1.5 while TiO2 can provide at 2.6. Refractive indexes in between are achievable with other materials such as SiN, HfO2 etc. All other parameters were kept the same as in the platelet radius sweep, i.e. the number of top and bottom DBR mirror pairs, the simulation box radius, and the PML thickness. The mesh was again a triangular- quadrilateral hybrid with maximum element size h ∼ 1/12×λref/nlocal. The results from the simulations are presented in Figures 4.5, a zoomed in version in 4.6 and a 3D representation in 4.7. 1.6 1.8 2 2.2 2.4 2.6 0 1000 2000 3000 4000 5000 6000 Figure 4.5: Threshold material gain as a function refractive index and thickness of the cladding. 35 4. Results 2 2.1 2.2 2.3 2.4 2.5 2.6 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Figure 4.6: Zoomed in image at the local minima in threshold material gain versus refractive index of the cladding. 300 250 200 0 1.6 150 1.8 2000 2 2.2 1002.4 4000 2.6 6000 8000 Figure 4.7: 3D representation of the choice of cladding. 36 4. Results The simulation results imply the existence of a global minimum around 2.3 < nClad < 2.5. For nClad close to global minimum, the threshold is lowest for a thick cladding layer. For 1.5 < nClad < 2.2, the cladding thickness has a large effect on the threshold as it seem to vary in a range between 2000 cm−1 < gth < 4000 cm−1. From the 2D perspective (Figure 4.5), the variations seem random, but in 3D space (Figure 4.7), there are clearly some smooth fluctuations of threshold in the nClad − RClad- plane. Furthermore, focusing on the global minimum, a thicker cladding layer seem to lower the threshold. This trend saturates at RClad = 260 nm, thus further thickening of the cladding layer won’t affect the threshold much. The lowest threshold is gth ∼ 500 cm−1, which was obtained for nClad = 2.44 and a thickness of RClad = 280 nm. A natural question that comes into mind is to what extent these optimal pa- rameters can be realized. A thickness of a few 100 nm is possible to deposit with current techniques. Undoped GaN has a refractive index of around 2.5 and could potentially be grown as a shell structure around the platelet. Another possibility is to use dielectric TiO2 that has a refractive index of around ∼ 2.4− 2.5. The choice of cladding material, however, affects the complexity in the fabrication, which must be taken in consideration when designing the cladding layer. 4.3 Concave DBR structure The wide angle scattering losses deriving from the pitched top section of the platelet is one of the main causes of the high threshold material gain. It is therefore necessary to find a way to compensate for those, by reducing the scattering angle. A concave mirror has the ability to focus a divergent light beam. Thus, a concave DBR might be able to reduce the threshold. A platelet radius sweep was performed for the concave top DBR structure de- scribed in Figure 4.8. To model a concave top DBR structure is challenging. Firstly, the concave structure complicates the generation of the hybrid mesh. This was solved Figure 4.8: Comparison between flat top DBR and a concave top DBR. 37 4. Results by reducing the maximum element size by 10 % to h ∼ (1 − 0.1)1/12 × λref/nlocal. Secondly, the concave structures spread out in the radial direction, thus they end up intersecting with the PML layers. To avoid oblique interfaces at the PML bound- ary, the simulation box size was adjusted so that all bent structures fit inside the simulation box. The simulation box sizes used for the platelet radius simulations are presented in Figure 4.9. 200 400 600 800 1000 1200 1400 1600 2000 2200 2400 2600 2800 3000 3200 3400 3600 Figure 4.9: Adjusted simulation box sizes for concave top DBR structures, to avoid oblique interfaces at the PML boundary. The simulation box sizes are rounded to closest multiple of 100 nm. Except for the maximum element size and the simulation box size, the other parameters were kept consistent with the previous platelet radius sweeps described earlier. The results from the platelet radius sweep of the concave top DBR structure is compared with the previous sweeps with flat DBRs in Figure 4.10. 900 1000 1100 1200 1300 1400 1500 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Figure 4.10: Comparison between flat and concave top DBR structures. 38 4. Results A concave top DBR leads to reduced losses and thus a significantly lower thresh- old material gain compared to having a flat top DBR on a pitched platelet. For in- stance atRpl = 1300 nm the threshold obtained with concave top DBR is 1713.9 cm−1 which is almost half to that of 3145.7 cm−1 for the flat top DBR case. Also, note that the threshold values for the concave top DBR are not seldom even comparable to the straight top platelet case, i.e. below 2000 cm−1. However, as seen in Figure 4.10 the threshold does not vary smoothly between different platelet radius points. In order to conclude if these fluctuations are real or related to numerical artifacts, an additional sweep has been performed with a smaller step size of 2.5 nm (compared to 25 nm in Figure 4.10). The results are presented in Figure 4.11. As seen in the figure, the finer sampling shows a continuous and smooth relation between threshold and platelet radius. Consequently, it is concluded that the fluctuations do not come from numerical noise and are actually part of the trend. 1250 1260 1270 1280 1290 1300 1310 1320 1330 1340 1350 1360 1370 1380 0 1000 2000 3000 4000 5000 6000 Figure 4.11: Comparison between flat and concave top DBR structures. 4.3.1 Cladding thickness and refractive index To further study the effects of the concave top DBR, a cladding sweep similar to the one presented in Section 4.2 was performed. In order to compare these results with the results from the flat top DBR case (Figure 4.5), the platelet radius was again pinned at Rpl ∼ 1300 nm. Except for the maximum element size and the simulation box size, the other parameters were kept consistent with the previous cladding sweeps described earlier. The results from the simulation are presented in Figure 4.12 39 4. Results 1.6 1.8 2 2.2 2.4 2.6 0 1000 2000 3000 4000 5000 6000 7000 8000 Figure 4.12: Threshold material gain as a function of refractive index and thickness of the cladding, for a concave top DBR structure. As implied by the figure, the cladding effect on the threshold material gain for a concave top DBR case is completely different from the flat top DBR case in Figure 4.5. The trend in Figure 4.12 has no sweet-spot valley and is instead monotonically increasing for higher cladding refractive index. Hence, in the case of a concave top DBR, it seems favorable to have a cladding with low refractive index. However, the figure also shows a sharp dip in threshold at around 2.5 for thick cladding, at similar threshold levels as for low refractive indices. Due to coarse sampling there might be some interesting features hidden between the 2.4 and 2.5 sampling points. Due to limited time, a finer sampled cladding sweep is left as future work, proceeding this thesis. 4.4 Number of DBR mirror pairs In all previous simulations, the number of top DBR mirror pairs was 15. To check the influence of mirror pairs on the threshold, this number was varied between 12 up to 25. The results are presented in Figure 4.13 and shows that the threshold can be further reduced by adding more mirror pairs. However, 18 mirror pairs is sufficient, hence adding more mirror pairs won’t affect the threshold by any significant amount. By changing the number of mirror pairs from 15 to 18, the threshold material gain for the concave case dropped from 1714 cm−1 to 1682 cm−1, while the flat case dropped from 3116 cm−1 to 3057 cm−1. 40 4. Results 12 13 14 15 16 17 18 19 20 21 22 23 24 25 1500 2000 2500 3000 3500 4000 Figure 4.13: Threshold material gain as a function of the number of mirror pairs in a flat and a concave top DBR. 4.5 Active region alignment The vertical alignment between the active region and the field anti-node is slightly off in some cases (see Figure 4.14 as typical example). The significance of miss- alignment was investigated by moving the active region slightly. In the simulations both the MQW and the EBL were displaced by δactive with δactive > 0 referring to positive z-direction. The case studied had a concave top DBR with 18 mirror pairs to further reduce the threshold material gain. The platelet radius was Rpl = 1300 nm, and the cladding thickness was RClad = 100 nm with a refractive index of nClad = 1.5 . The results are presented in Figure 4.15. The figure shows that a minimum threshold material gain is achieved if the active region is actively displaced by δactive = −4 nm from its original position, which confirms the miss-alignment implied in Figure 4.14. By looking at the trend, the minimum region is fairly wide. Hence, small displacements up to ±5 nm from the threshold minimum have insignificant effect on the threshold material gain. Thus the design is robust in regard to small variations in active region placement. For larger displacements the overlap between the active region and the field anti-node becomes poor, and the resulting threshold material gain undoubtedly increases. 41 4. Results Figure 4.14: Vertical cross-section that showcases the electric field alignment with active region. -20 -16 -12 -8 -4 0 +4 +8 +12 +16 +20 1500 2000 2500 3000 3500 4000 Figure 4.15: Threshold material gain as a function of active region displacement from default position. 42 4. Results 4.6 Absorption effects The influence of absorption was briefly studied for a case that from the above men- tioned simulations exhibits a low threshold. A platelet-based laser structure with a concave top DBR of 18 mirror pairs; a platelet radius of Rpl = 1300 nm, and a 100 nm thick cladding with refractive index 1.5, was chosen as candidate. The active region was shifted by δactive = −4 nm for improved alignment. Absorptive properties were added to the indium tin oxide (ITO) and to the GaN through the extinction coefficient. The absorption coefficient for ITO was set to αITO = 1000 cm−1 and to αGaN = 10 cm−1 for GaN. The resulting threshold material gain was Without absorption: gth = 1624 cm−1 With absorption: gth = 2465 cm−1 (4.1) hence, absorption will significantly increase the threshold material gain. 4.7 Reliability and error estimates Early results1 based on the beam propagation method BPM showed implications that platelet-based GaN VCSELs with small radii exhibit extremely high thresholds. From the platelet radius simulations presented in this thesis it is implied that larger radii is necessary to achieve reasonable threshold levels. In order to validate these results, one BPM simulation was performed for Rpl = 1300 nm with flat top DBR. The output threshold material gain from the BPM simulation was 2600 cm−1, which is in fairly good agreement with the FEFD simulation results that gave 3100 cm−1. In the BPM simulation, the bottom part of the platelet had rounded corners, unlike the platelet model in the FEFD simulation. These rounded corners suppos- edly have similar focusing ability as the implementation of a concave top DBR has, thus a significant reduction of threshold material gain is expected. Hence, this might explain the difference in threshold material gain results above. For a more rigor- ous validation of the results from the FEFD simulations, a BPM simulation with straight corners should have been used, but due to lack of time this has been left for future work. However, since both the BPM and FEFD gave similar results, one can still say with confidence that the reliability of the gain sweep method has been encouraged. In order to test the accuracy of the results presented, an error estimate has been made. All error estimate simulations have been performed on a representative structure with Rpl = 1300 nm, Rclad = 100 nm and nclad = 1.5. Both flat and concave top DBRs have been considered. The number of mirror pairs were 15 for the top DBR and 10 for the bottom DBR, unless otherwise stated. Only a realistic tilt angle of θ = 58 ◦ was considered, and the active region was left in default position. The simulation box size was per default 2000 nm for the flat top DBR case and 3200 nm for the concave top DBR case. 1Photonics Laboratory, Department of Microtechnology and Nanoscience, Chalmers University of Technology 43 4. Results 4.7.1 Simulation Box Radius The simulation box radius was swept from 2000 nm up to 13000 nm for both a flat and a concave top DBR case. The simulation results are shown in Figure 4.16. The simulation box is expected to have insignificant effect on the threshold material gain, thus a mean and a 2σ-standard deviation have been calculated. The threshold material gain varies within less than 1 % from the mean value. Hence, the results presented in this thesis are robust with respect to simulation box size. Concave DBR: gth = 1703, 6± 5, 5 (2σ, 0, 32 %) [1/cm] Flat DBR: gth = 3106, 5± 20, 6 (2σ, 0, 66 %) [1/cm] (4.2) 0 3000 6000 9000 12000 15000 1680 1685 1690 1695 1700 1705 1710 1715 1720 0 3000 6000 9000 12000 15000 3060 3070 3080 3090 3100 3110 3120 3130 3140 Figure 4.16: Threshold material gain dependency on simulation box radius. 4.7.2 PML Thickness The thickness of the perfectly matched layer (PML) was swept from 200 nm up to 420 nm for both a flat and a concave top DBR case. The simulation results are shown in Figure 4.17. If the PML is tuned properly, it is expected to have insignificant effect on the threshold material gain. However, the simulation results imply that the threshold material gain first decreases with thicker PML, and then flattens out. These results suggest that a thicker PML than the one used in the previous simulations would have been preferred for more robust results. Howbeit, the effect of implementing a thicker PML is minuscule as the threshold material gain would only change by less than 1 %. Hence, the results presented in this thesis can also be treated as robust with respect to PML thickness. 44 4. Results Furthermore, the mean values presented in Equation (4.3) are somewhat lower than the threshold material gains presented in Equation (4.2). The bias is due to the top DBR had 18 mirror pairs instead of 15 as default. Concave DBR: gth = 1680, 5± 4, 7 (2σ, 0, 28 %) [1/cm] Flat DBR: gth = 3050, 4± 11, 4 (2σ, 0, 36 %) [1/cm] (4.3) 200 250 300 350 400 1660 1665 1670 1675 1680 1685 1690 1695 1700 200 250 300 350 400 3000 3010 3020 3030 3040 3050 3060 3070 3080 Figure 4.17: Threshold material gain dependency on PML thickness. 100 % 90 % 80 % 70 % 60 % 1680 1685 1690 1695 1700 1705 1710 1715 1720 100 % 90 % 80 % 70 % 60 % 3060 3070 3080 3090 3100 3110 3120 3130 3140 Figure 4.18: Threshold material gain dependency on maximum element size. 45 4. Results 4.7.3 Maximum Element Size The maximum element size was swept from h ∼ λref/12nlocal(100 %) down to 60 % for both a flat and a concave top DBR case. The simulation results are shown in Figure 4.18. If the mesh is good enough, a reduction of maximum allowed element size would neither affect the field solution nor the output threshold material gain. As can be seen in the figure, the threshold material gain is very robust with respect to the maximum element sizes simulated (2σ < 0.1 %), hence a maximum element size of h ∼ λref/12nlocal is sufficient to obtain reliable results. Concave DBR: gth = 1704, 2± 1, 4 (2σ, 0, 08 %) [1/cm] Flat DBR: gth = 3116, 37± 1, 6 (2σ, 0, 05 %) [1/cm] (4.4) 46 5 Conclusions The finite element frequency domain method has been used to simulate the optical losses in a novel nanowire-based GaN vertical-cavity surface-emitting laser. By using a gain sweep method, the threshold material gain has been calculated for a number of different design parameters. The radius of the platelet is the most crucial parameter when it comes to diffraction losses. The simulations show that for small radii, comparable to the resonance wavelength of the cavity, diffraction losses are inevitable. Hence, in order to reach decent threshold material gain and achieve lasing with the current geometric design, the platelet radius must be increased to at least 1200 nm to 1300 nm. Furthermore, the pitch of the platelet introduces additional diffraction losses, which can be suppressed by depositing a conformal top DBR, i.e. a concave mirror with focusing features. For further reduction of threshold material gain, the number of mirror pairs in the top DBR can be increased. By increasing from 15 to 18, the threshold material gain is reduced by about 50 cm−1. Adding even more mirror pairs does not reduce the threshold significantly further. Furthermore, the active region is somewhat miss- aligned with the optical field anti-nodes in the original structural design. Optical absorption in the materials clearly has a significant effect on the thresh- old material gain, and for the particular case studied here adding realistic absorption in the ITO and GaN-layers increases the threshold from 1600 cm−1 to 2500 cm−1. Hence this effect should be investigate more rigorously in future work. The thickness and materials of the cladding layer also strongly affects the threshold material gain, but the influence of the cladding layer is clearly different between having a flat top DBR or a concave top DBR. The results from the flat top DBR case imply that a small index contrast between the core and cladding materials is beneficial, while the results from the concave top DBR case show a completely opposite trend. Thus, further studies of the influence of threshold material gain from the cladding layer properties are required. Error estimation simulations were performed to check the sensitivity of the method to simulation box size, PML thickness, and maximum element size in the mesh. When these three parameters were increased by four times, two times, and reduced by 60 %, respectively, the threshold material gain only changed with < 1 %, demonstrating the robustness of the method. To validate the method one simu- lation point was compared to beam propagation method (BPM) for the case of a large platelet with a radius of 1300 nm. The BPM gave a threshold material gain of 2600 cm−1, which is comparable to 3100 cm−1 given by the FEFD method. 47 5. Conclusions To summarize, an FEFD method had been used to calculate the threshold ma- terial gain in novel nanowire-based GaN vertical-cavity surface-emitting lasers. The most critical parameter in the design is the platelet radius, where it is shown that the radius should be at least 1200 to 1300 nm to achieve low enough threshold ma- terial gain. The cladding also has a strong effect on the threshold material gain, but further investigations are required to understand its complex influence. The robustness of the method and comparable results to a BPM method validates its accuracy. Being far less computationally heavy than BPM, it will be an important tool in the future design of novel lasers. 48 Bibliography [1] Å. Haglund, E. Hashemi, J. Bengtsson, J. Gustavsson, S. Carlsson, M. Stattin, M. Calciati, and M. Goano: Progress and challenges in electrically pumped GaN-based VCSELs. Proc. of SPIE, 9892, 2016 [2] S. Arafin, X. Liu, and Z. Mi: Review of recent progress of III-nitride nanowire lasers. Journal of Nanophotonics, 7, 2013 [3] A. Einstein: Zur Quantntheorie der Strahlung. Physikalische Zeitschrift, 18, pp. 121-128, Mar 1917 [4] K. Thyagarajan, and A.K. Ghatak: Lasers Theory and Applications. Plenum Press, 1981 [5] S. M. Sze, and KWOK K. Ng: Physics of Semiconductor Devices. John Wiley & Sons, Inc., Third Edition, 2007 [6] G. K. Svedensen, H. Weman, and J. Skaar: Investigations of Bragg reflectors in nanowire lasers. Journal of Applied Physics, 111, 2012 [7] R.