Numerical Self-Consistent Analysis of VCSELs

Abstract

Vertical-cavity surface-emitting lasers (VCSELs) yield single-longitudinal-mode operation, low-divergence circular output beam, and low threshold current. This paper gives an overview on theoretical, self-consistent modelling of physical phenomena occurring in a VCSEL. The model has been experimentally confirmed. We present versatile numerical methods for nitride, arsenide, and phosphide VCSELs emitting light at wavelengths varying from violet to near infrared. We also discuss different designs with respect to optical confinement: gain guidance using tunnel junctions and index guidance using oxide confinement or photonic crystal and we focus on the problem of single-transverse-mode operation.

1. Introduction

Currently there are two distinctly different classes of Fabry-Perot semiconductor diode lasers: edge-emitting lasers (EELs) and vertical-cavity surface-emitting lasers (VCSELs). Because of the details of their structure, VCSELs have a number of unique features that distinguish them from conventional EELs [1]: inherent single-longitudinal-mode operation, low-divergence nonastigmatic circular output beams, low threshold current at room-temperature (RT) continuous-wave (CW) operation, device geometry suitable for integration into two-dimensional laser arrays, compatibility with vertical-stacking architecture, the ability to be modulated at very high frequencies, and the possibility of in situ testing.

While EELs usually emit many longitudinal modes around the maximal optical gain wavelength, VCSELs emit a single longitudinal mode at the wavelength determined by the cavity design. Therefore, EEL cavities are always tuned to their maximal optical gain but those of VCSELs may be intentionally detuned, which gives an additional degree of freedom to the VCSEL design. As a result, designers of EELs can propose devices emitting radiation of wavelength solely determined by their activeregions. In the case of VCSELs, however, it is possible to design a device emitting radiation at wavelength somewhat different (practically always longer) from that associated with their active-region structure.

Computer simulations of laser operation enable anticipating the laser performance in more efficient and inexpensive way than the trial-and-error method. However, VCSEL modelling is a very involved task because of its multilayered structure (sometimes containing as many as several hundred layers) often of nonplanar or buried-type architecture, with many heterojunctions, graded layers, strained layers, quantum wells (QWs), quantum dots (QDs) or quantum wires, superlattices, oxide and oxidized layers, or mesa structures, and so forth. Therefore, in advanced VCSEL models, all material and structure parameters should be functions of the local material composition. Additionally, physical phenomena taking place during the VCSEL operation, that is, optical, electrical, thermal, and mechanical phenomena, are interrelated. These interactions create a network of mutual interrelations. Therefore comprehensive VCSEL modelling should have self-consistency, where as many as possible interactions between individual physical phenomena taking place within the VCSEL should be accounted for.

Following the general principles given in [2], the comprehensive self-consistent VCSEL model should contain the following mutually interrelated parts.(1)A three-dimensional (3D) optical model describing, for successive cavity modes, their modal gain and loss, lasing thresholds, emission wavelengths, and optical field distributions within the laser cavity.(2)A 3D electrical model characterizing the current spreading between the top and the bottom contacts through the centrally located active region, the injection of carriers of both polarity into the QW or QD active region, and their subsequent radiative or nonradiative recombination after radial out-diffusion within the active region and the possible overbarrier carrier leakage.(3)A thermal model characterizing generation of a heat flux (nonradiative recombination, reabsorption of spontaneous radiation as well as volume and barrier Joule heating), its 3D flow from heat sources towards a heat sink, and its 3D spreading within a heat sink.(4)A recombination model describing phenomena within the QW or QD active region, that is, furnishing information about the optical gain process being the result of radiative bimolecular recombination (to enable determination of the active-region optical gain spectrum) as well as nonradiative monomolecular and Auger recombinations.(5)a piezoelectric model (important in nitride materials) first describing mechanical stresses within the VCSEL structure necessary to determine the piezoelectric polarization and then introducing piezoelectric and spontaneous polarization effects leading to some modifications of potential distributions especially important in quantum wells (QWs).

Additionally in a comprehensive VCSEL model, as many as possible linear and nonlinear interactions between the above individual optical, electrical, thermal, recombination, mechanical, and piezoelectric phenomena should be taken into account with the aid of the self-consistent approach including the effects of(i)thermal focusing, that is, the temperature dependence of refractive indices,(ii)self-focusing, that is, dependence of refractive indices on carrier concentration,(iii)gain-induced wave guiding, that is, the temperature, carrier concentration, and wavelength dependences of the extinction coefficient,(iv)temperature (and sometimes also carrier concentration) dependence of thermal conductivities,(v)temperature and carrier concentration dependences of electrical conductivities,(vi)the temperature, carrier concentration, and wavelength dependences of optical gain,(vii)temperature and carrier concentration dependences of the energy gaps,(viii)quantum-confined Stark effect (QCSE), that is, an impact of spontaneous and piezoelectric polarization on the QW potential distribution leading to reduced overlapping of the wave functions of electrons and holes trapped within the QW, raising both their QW levels and carriers enhanced escape possibilities.

Truly comprehensive VCSEL models constitute an advanced and integrated tool for studying all aspects of VCSEL performance in the whole complexity of all interrelated physical phenomena taking part in their operation. 3D profiles of all model parameters within the whole device volume should be determined in each calculation loop not only by the various chemical composition of the layers, but by the self-consistent calculation algorithm which takes into account current 3D profiles of temperature, current density, carrier concentration, mode radiation intensity, and mechanical stresses. Reaching such self-consistency is especially important for modelling high-power and/or high-temperature VCSEL operation. Such a comprehensive VCSEL model has been reported by Sarzała and Nakwaski [3] and Sarzała [4], and its validity has been confirmed experimentally for QD VCSELs in [5]. Above the lasing threshold, the spatial-hole-burning (SHB) effect should also be taken into account [6–8]. A local increase in a lasing intensity of the already excited lowest-threshold mode enhances (as a result of a stimulated emission) recombination of active-region carriers leading to a decrease in their concentration and a decrease in the local optical gain at places of high mode radiation intensity. The modal optical gain of this lowest-threshold mode saturates (or may even be reduced) which may favour excitation of the next transverse mode of an intensity spatial profile better suited to this new optical gain profile. Therefore, the above-threshold analysis of a VCSEL operation requires considering the SHB effect.

Similar to a chain being as strong as its weakest link, a model is as exact as its least precise part. Therefore, it is useless to improve the more accurate parts of the model and not to care about less exact portions. All parts of the model should exhibit similar accuracy. In any individual case, a reasonable compromise should be reached between high modelling fidelity and its practical convenience depending on a main modelling goal, importance and urgency of expected results, available equipment, and also financial possibilities. Rigorous modelling from the first principles is usually very time-consuming and often requires powerful super-computers equipped with enormous memory. Sometimes such modelling is justified but very often more approximate (but quicker and cheaper) models are acceptable and provide important technological suggestions. The computational cost of applying a scalar or a vectorial optical model, a simplified or a more involved drift-diffusion electrical model, Fermi’s Golden Rule or a many-body optical gain determination, and so forth, has to be compared before the modelling approach is chosen. Additionally, it should be taken into consideration that rigorous models are very time-consuming, which practically excludes their application in self-consistent modelling. Often, modelling errors due to the lack of self-consistency are worse than those created by using a simpler approximate method.

The main goal of the current paper is to present versatile numerical methods for nitride, arsenide, and phosphide VCSELs emitting light at wavelengths varying from violet to near infrared. We also discuss different device designs with respect to optical confinements: gain guidance using a tunnel junction and index guidance using an oxide or Photonic Crystal. We also restrict our discussion to the problem of enhanced single transverse-mode operation. Such operation is required in several VCSEL applications as, for example, telecommunications and gas detection. Single mode operation can only be achieved for small aperture in standard index-guided VCSELs, which limits the available single-mode output [5]. Higher-output large-size index guided VCSELs usually exhibit multimode operation, especially at higher temperatures. This results from the considerable current-crowding effect near the active-region edges, which—despite the smooth radial carrier diffusion profile in the active region—still favours higher-order transverse modes [5, 8]. In gain-guided VCSELs, on the other hand, single-mode operation is possible for broader active regions, at the expense of higher threshold currents. So as one can see, the requirements of higher-power single-mode operation and low lasing threshold usually contradict each other. Nevertheless, various approaches of enhancing the single-mode VCSEL output power have been reported [8–14].

Even if the VCSEL is single mode, there still can be observed polarization switching between two orthogonal polarisations of the fundamental mode. In a number of applications such polarization switching cannot be tolerated. For this reason, the polarization stabilization is an important issue in the design of modern VCSELs. Traditional, commonly used VCSELs have geometries with cylindrical symmetry and are usually grown on (100)-oriented substrates. Thus, they have no a priori defined polarization, that is, the LP-polarized modes can be oriented in any direction perpendicular to the VCSEL symmetry axis. In practice, there is always some small anisotropy induced by residual strain introduced during fabrication and the electrooptic effect. Therefore, the electric field of the fundamental mode (LP mode) is almost always linearly polarized along [011] or [011] crystallographic axis, while the one of the first-order modes is perpendicular [15]. During the operation of the laser, a change of temperature, current, or stress can induce polarization switching, that is, the VCSEL starts emitting light in the polarization orthogonal to the previous one. This effect can be a subject to bistability, hysteresis or random dynamic behavior [16, 17]. Furthermore, optical feedback [18] or optical injection [19] can introduce even more polarization instabilities in the system. There exist several methods for providing polarization control in VCSELs. They can be attributed to one of the following groups: introduction of anisotropic gain, use of noncircular resonators or polarization-dependent mirrors, an external feedback, and finally application of photonic crystals. Here, we are focused on the numerical design of the optimal photonic crystal aimed for stabilization of polarization in VCSELs.

The paper is organized as follows. Section 2 describes the theoretical self-consistent model of physical phenomena taking part in VCSEL, and Section 3 concerns the experimental validation of the model based on a 1300 nm QD VCSEL. Comparison of scalar and vectorial optical models for 400 nm nitride VCSEL is presented in Section 4. Optimization of oxide-confined VCSELs is presented in Section 5. Design rules of photonic crystal VCSELs for enhanced single-mode and single-polarization operation are given in Sections 6 and 7, respectively. Finally, in Section 8 we conclude.

2. Numerical Model

The simulated performance characteristics of VCSELs are determined with the aid of numerical models presented in this section. Typically, personal computers are used in simulation process, and the resulting calculation time to achieve a single result of any characteristic varies from several minutes to one day. The shorter calculation times arise from using 1D or 2D electrical and thermal models combined with 2D scalar optical approach, while the longer duration occurs when using all 3D models and vectorial optics.

2.1. Optical Models

2.1.1. Scalar Optics

The scalar optical approach for VCSEL simulation was introduced by Hadley [20] and latter modified by Wenzel and Wünsche [21]. In this effective-frequency method (EFM) the optical field 𝐸(𝑟,𝑧,𝜑) is assumed to be linearly polarized and spatially separated:𝐸𝑟,𝑧,𝜑=𝐸_𝑟,𝑧(𝑟,𝑧)𝐸_𝑟(𝑟)exp𝑖𝐿𝜑, (1) where 𝑧 is the direction of the light propagation in the VCSEL, 𝑟 is the radial coordinate in the plane of epitaxial layers, φ is the azimuthal coordinate, and 𝐿=0,1,2,… is the azimuthal mode number. The above enables reducing the wave equation within the VCSEL resonator to the following two mutually interrelated nearly-one-dimensional wave equations along the axial and the radial VCSEL directions:𝑑²𝑑𝑧²+𝑘²₀𝜀𝐸_𝑟,𝑧=𝑣_eff(𝑟)𝑘²₀√𝜀𝑛_𝑔𝐸_𝑟,𝑧,𝑑²𝑑𝑟²+1𝑟𝑑−𝐿𝑑𝑟𝑟²²+𝜈_eff(𝑟)𝑘²₀√𝜀𝑛_𝑔𝑟𝐸^𝐿_𝑟(𝑟)=𝜈𝑘²₀√𝜀𝑛_𝑔𝑟𝐸^𝐿_𝑟(𝑟). (2) For simplicity, we assigned 𝜀≡𝜀(𝑟,𝑧) and 𝐸_𝑟,𝑧≡𝐸_𝑟,𝑧(𝑟,𝑧), 𝜈_eff to represent the effective frequency, 𝑘₀ is the vacuum wave number, ε is the dielectric constant of the layer, and 𝑛_𝑔 is the group index. The boundary conditions (BCs) are realized by the introduction of the region in which the solutions of the (2) are assumed as decaying functions [22]. The last layer that is responsible for decaying of the field in the radial direction is placed 20 μm from the edge of the active region. The dimensionless complex parameter 𝜈 plays the role of an eigenvalue.

2.1.2. Vectorial Optics

The computational efficient vectorial approach to VCSELs described in this section is based on the plane wave admittance method (PWAM) [19]. The main objective of the method is the transformation of Maxwell’s equations to a form of the characteristic equation, which will be utilized in the further vectorial analysis of the VCSEL. To this aim, the Cartesian coordinate system is oriented in such a way that the x-y plane is parallel to the epitaxial layers, whereas the 𝑧 direction is the direction of the light propagation. Then, Maxwell’s equations can be expressed in the form∇×𝐄𝑥,𝑦,𝑧,𝑡=−𝜇𝜇₀𝜕𝐇𝑥,𝑦,𝑧,𝑡𝜕𝑡∇×𝐇𝑥,𝑦,𝑧,𝑡=𝜀𝜀₀𝜕𝐄𝑥,𝑦,𝑧,𝑡𝜕𝑡 (3) with μ, 𝜇₀ being the magnetic permitivity diagonal tensors for the material and vacuum, respectively, and ε, 𝜀₀ being the respective diagonal tensors of dielectric constant. E and H are the vectors of the electric and magnetic fields. We further assume that(1)harmonic time dependence of the fields is of the form ∼exp(𝑖𝜔𝑡), with ω being the angular frequency of the propagating wave in the vacuum,(2)the structure consists of uniform (in the propagation z-direction) parallel layers, which yields Maxwell’s set of the equations in the static form∇×𝐄𝑥,𝑦,𝑧=−𝑖𝜔𝜇𝜇₀𝐇𝑥,𝑦,𝑧,∇×𝐇𝑥,𝑦,𝑧=𝑖𝜔𝜀𝜀₀𝐄𝑥,𝑦,𝑧. (4) Eliminating the z-components of the electric and magnetic fields from the above equations results in𝜕²_𝑧𝐸_𝑥𝐸_𝑦=1𝜔²𝜀₀𝜇₀−𝜕_𝑥1𝜀_𝑧𝜕_𝑥−𝜔²𝜀₀𝜇_𝑦𝜇₀𝜕_𝑥1𝜀_𝑧𝜕_𝑦−𝜕_𝑦1𝜀_𝑧𝜕_𝑥𝜕_𝑦1𝜀_𝑧𝜕_𝑦+𝜔²𝜀₀𝜇_𝑥𝜇₀×−𝜕_𝑦1𝜇_𝑧𝜕_𝑦−𝜔²𝜇₀𝜀_𝑥𝜀₀𝜕_𝑦1𝜇_𝑧𝜕_𝑥−𝜕_𝑥1𝜇_𝑧𝜕_𝑦𝜕_𝑥1𝜇_𝑧𝜕_𝑥+𝜔²𝜇₀𝜀_𝑦𝜀₀𝐸_𝑥𝐸_𝑦. (5) The electromagnetic fields as well as the magnetic and electrical permitivities are decomposed in orthonormal complete basis of exponential functions:Φ_𝑢=^𝑁_𝑛,𝑚Φ_𝑢^𝑛,𝑚𝜑_𝑛,𝑚,𝜂_𝑢=^∞_𝑛,𝑚𝜂_𝑢^𝑛,𝑚𝜑_𝑛,𝑚, (6) where Φ_𝑢 are the arbitrary field components of the electric or magnetic field, 𝜂_𝑢 are the components of magnetic or electric permitivity, and 𝑢=𝑥,𝑦,𝑧. The basis functions have been defined in the form of product of two functions, which satisfy the orthonormality and completeness of the basis:𝜑_𝑛,𝑚=exp𝑖2𝜋𝑛𝐿_𝑥+𝑘_𝑥𝑥+𝑖2𝜋𝑚𝐿_𝑦+𝑘_𝑦𝑦, (7) where 𝐿_𝑥 and 𝐿_𝑦 correspond to the dimensions of the calculation window along the 𝑥 and 𝑦 axis and 𝑘_𝑥 and 𝑘_𝑦 are corresponding components of wavevector in the x-y plane. Using these assumptions modifies the set of (5) to the form in which fields and permitivities are replaced with coefficients of exponential expansions:𝜕²_𝑧𝐸^𝑛_𝑦𝐸^𝑛_𝑥1=−𝑘²₀(𝐆+𝐤)𝐆^′′+𝐤𝜅^𝐆′′_𝑧^−𝐆−𝑘²₀𝜇^𝐆′′_𝑥^−𝐆(𝐆+𝐤)𝐆^′′+𝐤𝜅^𝐆′′_𝑧^−𝐆−(𝐆+𝐤)𝐆^′′+𝐤𝜅^𝐆′′_𝑧^−𝐆𝑘²₀𝜇^𝐆′′_𝑦^−𝐆−(𝐆+𝐤)𝐆^′′+𝐤𝜅^𝐆′′_𝑧^−𝐆×(𝐆+𝐤)𝐆^′′+𝐤𝛾^𝐆′′_𝑧^−𝐆−𝑘²₀𝜀̃^𝐆′′_𝑦^−𝐆(𝐆+𝐤)𝐆^′′+𝐤𝛾^𝐆′′_𝑧^−𝐆−(𝐆+𝐤)𝐆^′′+𝐤𝛾^𝐆′′_𝑧^−𝐆𝑘²₀𝜀̃^𝐆′′_𝑥^−𝐆−(𝐆+𝐤)𝐆^′′+𝐤𝛾^𝐆′′_𝑧^−𝐆𝐸^𝐆_𝑦𝐸^𝐆_𝑥. (8) The boundary conditions assumed for the set of equations are fulfilled by absorbing perfectly matched layers [23]. The final set (8) can be written as𝜕²_𝑧𝐄=𝐑_𝐇𝐑_𝐄𝐄≡𝐐_𝐄𝐄. (9) In the same manner one arrives at the equation for magnetic field:𝜕²_𝑧𝐇=𝐑_𝐄𝐑_𝐇𝐇≡𝐐_𝐇𝐇. (10) Equations (9) and (10) can be solved in a base, which reduces the matrices 𝐐_𝐄 and 𝐐_𝐇 to diagonal forms:𝜕²_𝑧𝐄+𝚪²_𝐄𝜕𝐄=0,_2𝑧𝐇+𝚪²_𝐇𝐇=0, (11) where 𝐄 and 𝐇 stand for electric and magnetic fields in the new base and can be defined as𝐄=𝐓_𝐄⁻¹𝐄,𝐇=𝐓_𝐇⁻¹𝐇, (12) where the matrices 𝐓_𝐄 and 𝐓_𝐇 diagonalize 𝐐_𝐄 and 𝐐_𝐇, respectively:𝐓_𝐄^−𝟏𝐐_𝐄𝐓_𝐄=𝚪²_𝐄,𝐓_𝐇^−𝟏𝐐_𝐇𝐓_𝐇=𝚪²_𝐇. (13) The solution of (11) has the well-known form of a standing wave:𝐄(𝑧)=cosh𝑖𝚪_𝐄𝑧𝐀_𝐄+sinh𝑖𝚪_𝐄𝑧𝐁_𝐄,𝐇(𝑧)=cosh𝑖𝚪_𝐇𝑧𝐀_𝐇+sinh𝑖𝚪_𝐇𝑧𝐁_𝐇. (14) Additionally, the fields 𝐄(𝑧) and 𝐇(𝑧) and their z-derivatives satisfy the continuity condition on the interface between the layers, that is,𝐇₀^(𝑖)−𝐇_𝑑^(𝑖)=𝐲_𝟏^(𝑖)𝐲_𝟐^(𝑖)𝐲_𝟐^(𝑖)𝐲_𝟏^(𝑖)𝐄₀^(𝑖)𝐄_𝑑^(𝑖), (15) where subscripts 0 and 𝑑 correspond to the bottom and the top of the layer, respectively,𝐲_𝟏=𝜶_𝐇tanh⁻¹𝑖𝚪_𝐄𝑑,𝐲_𝟐=−𝜶_𝐇sinh⁻¹𝑖𝚪_𝐄𝑑,𝜶_𝐇=𝐓_𝐄⁻¹𝐑_𝐇𝐓_𝐇⁻¹𝚪_𝐄. (16) Transformation of the fields from layer 𝑖 to layer 𝑖+1 is given by𝐄₀^(𝑖+1)=𝐓_𝐄^(𝑖+1)−1𝐓_𝐄^(𝑖)𝐄_𝑑^(𝑖)≡𝐭_𝐄^(𝑖)𝐄_𝑑^(𝑖),𝐇₀^(𝑖+1)=𝐓_𝐇^(𝑖+1)−1𝐓_𝐇^(𝑖)𝐇_𝑑^(𝑖)≡𝐭_𝐇^(𝑖)𝐇_𝑑^(𝑖). (17) The above relations of the field transformation between and within the layers can be employed to determine the characteristic equation by forcing boundary conditions, that is, zeroing of the field on the upper-most and lower-most layer boundaries. This allows finding the characteristic values of the problem, that is, the effective emission wavelength of the VCSEL, as well as the characteristic vectors, which determine the distribution of the electromagnetic field within the structure. Furthermore, one can choose a matching interface, which should be placed at a maximum of standing wave to reduce numerical error. That leads to the following relation between electric fields in neighboring layers:𝐓_𝐇^(𝑚)𝐘^(𝑚)_up𝐄_𝑑^(𝑚)=𝐓_𝐇^(𝑚)𝐘^(𝑙)_down𝐄_𝑑^(𝑙), (18) where𝐘^(𝑖)=−𝐲_𝟐^(𝑖)𝐭_𝐇^(𝑖)𝐘^(𝑖−1)𝐭_𝐄^(𝑖)−1−𝐲_𝟏^(𝑖)−1𝐲_𝟐^(𝑖)+𝐲_𝟏^(𝑖), (19) that is, the matrices Y are determined from a recurrence relation starting from the bottom 𝐘^(𝑚)_up and the top 𝐘^(𝑚)_down of the VCSEL. Thus, (18) results in𝐓_𝐇^(𝑚)𝐘^(𝑚)_up𝐓_𝐄^(𝑚)−1−𝐓_𝐇^(𝑚)𝐘^(𝑙)_down𝐓_𝐄^(𝑙)−1𝐄=𝐘𝐄=0, (20) which is the characteristic value equation. The solution of (20) determines the complex emitted wavelength. The imaginary part of wavelength determines the modal gain of the propagating wave.

2.2. Electrical Model

We utilize a finite-element (FE) electrical model of the VCSEL, which takes into account in a natural way the interplay of drift and diffusion currents for both electrons and holes. Current spreading can be determined from a three-dimensional (3D) potential distribution 𝑉(𝑟,𝑧,𝜑) by solving the Laplace equation:div𝜎𝑟,𝑧,𝜑grad𝑉𝑟,𝑧,𝜑=0, (21) where σ denotes the 3D profile of electrical conductivity and 𝑟, 𝑧, and φ represent a cylindrical coordinate system with the 𝑧-axis directed along the device axis. From now on, the following notation 𝑟⃗=(𝑟,𝑧,𝜑) will be used. For all the layers of the laser structure (except the active region), the conductivity σ depends on the material composition and doping, as well as on local temperature and local carrier concentration:𝜎(𝑟)⃗=𝑒𝑛(𝑟)⃗𝜇(𝑟)⃗. (22) In (22), 𝑒 is the electron charge, 𝑛 the carrier concentration, and μ the carrier mobility.

Generation and recombination phenomena within the active region provide the nonzero right-hand side of (21) (known as the Poisson equation in this case). Their relative influence is difficult to analyze theoretically. Therefore, they will be taken symbolically into account with the aid of an effective conductivity 𝜎_pn of the active-region material. Its value can be found using the differential Ohm law and the classical diode equation:𝜎_pn(𝑟)=𝛽_pn𝑗_pn(𝑟)𝑑_𝐴,𝐸ln𝑗_pn(𝑟)/𝑗_𝑠+1, (23) where 𝑗_pn is the p-n junction current density, 𝑗_𝑠 the reverse p-n saturation current density, and 𝑑_𝐴,𝐸 the cumulative active-region thickness including not only the active layers but also the barrier layers between them.

To obtain the 3D potential profile for the whole laser structure, it should be matched (using a self-consistent approach) with the aid of boundary conditions at all boundaries between the layers. The 3D current density 𝑗(𝑟)⃗ can then be found within the whole device volume from the differential Ohm law:𝑗(𝑟)⃗=−𝜎(𝑟)⃗grad[𝑉(𝑟)⃗]. (24) Thereafter, the carrier density profile 𝑛_𝐴(𝑟) within the active layer may be determined from the diffusion equation below threshold:𝐷_𝐴𝜕²𝑛_𝐴(𝑟)𝜕𝑟²+1𝑟𝜕𝑛_𝐴(𝑟)𝜕𝑟−𝐴𝑛_𝐴(𝑟)−𝐵𝑛²_𝐴(𝑟)−𝐶𝑛³_𝐴(𝑟)+𝑗_pn(𝑟)𝑒𝑑_𝐴=0, (25) where 𝐷_𝐴 is the ambipolar diffusion coefficient, (𝐴) is the monomolecular (on point defects) nonradiative recombination constant, (𝐵) the bimolecular recombination constant, and (𝐶) the Auger recombination coefficient.

2.3. Thermal Model

We utilize a FE thermal VCSEL model, for which the same mesh as the one generated for the FE electrical calculations is applied. The heat conduction equationdiv𝜆_𝑇(𝑟)⃗grad(𝑇(𝑟)⃗)=−𝑔_𝑇(𝑟)⃗ (26) may be solved in this way for the whole VCSEL structure. In the above equation, 𝜆_𝑇 stands for the 3D thermal conductivity coefficient and 𝑔_𝑇 for the 3D volume density of heat sources.

Nonradiative recombination and reabsorption of spontaneous radiation within the active region together with the volume Joule heating in all layers and the barrier Joule heating in the contacts are usually the main heat sources in the laser. The 3D heat flux spreading in the heat sink is determined by assuming that its external walls are kept at ambient temperature. With the exception of some specific cases, the mesa sidewalls and top laser facet are assumed to be thermally isolated because of the negligible heat abstraction mechanisms.

2.4. Overthreshold Model

The active-region carrier concentration distribution 𝑛_𝐴 is described by a time-independent (steady state) diffusion equation with adiabatic boundary conditions in one or two dimensions, that is, the active region is assumed to be so thin that 𝑛_𝐴 does not depend on 𝑧:𝜅Δ𝑛_𝐴(𝑟)⃗−𝐿𝑇,𝑛_𝐴+𝑆(𝑟)⃗=0, (27) where 𝜅 is the diffusion constant, 𝐿 are the carrier losses, which depend on the carrier concentration and local temperature (see (28)), and 𝑆 describes the carrier source (see (29)). The losses in the threshold case consist of three recombination processes: monomolecular, bimolecular, and Auger:𝐿𝑇,𝑛_𝐴=𝐴𝑛_𝐴+𝐵𝑛_𝐴²+𝐶𝑛_𝐴³, (28) where the recombination coefficients 𝐴, 𝐵, and 𝐶 may depend on the temperature. The sources are just the carriers injected into the active region of thickness 𝑑:𝑆()=𝑗𝑟⃗^⟂𝑒𝑑, (29) where 𝑗_⟂ denotes the current density in the direction perpendicular to the active region.

Now we want to add to our threshold model the ability to describe phenomena connected with the presence of a strong electromagnetic field, which causes the stimulated recombination of the carriers in the active region (where the optical gain is positive). This phenomenon is called spatial hole burning. We will simulate it by adding another (negative if the gain is positive) carrier source in the carrier diffusion equation. Since our equation assumes a steady state, we must assume that the stimulated recombination rate does not depend on time. This requires that the modes are stable—they do not change in time. This assumption may not be valid in many edge emitting lasers, but in VCSELs it is very natural. In case of single-mode operation, the stimulated emission losses are described by the following formula:𝐿_st𝑛_𝐴𝑛,𝑇,𝑥,𝑦=𝑃2𝑔_𝐴,𝑇,𝑟⃗𝑀(𝑟)⃗(1−𝑅)ℏ𝜔∫_𝑆𝑀, (30) where 𝑃 is the output power, 𝑔 is the optical gain, 𝑀 the mode profile, 𝑅 is the reflectivity of the output mirror (the other one is assumed to reflect 100%), and ℏω is the photon energy. Adding the above term to (27) results in the above-threshold carrier diffusion equation. For a given voltage 𝑈, which determines the temperature and carrier distribution, we find such a value of 𝑃, which gives such a carrier (and hence optical gain) distribution that the total modal gain is equal to 0.

In case of the multimode operation, stimulated emission losses for each mode are added. Each mode has its own power, gain, mode distribution, and photon energy. In this case we formally need to solve a system of 𝑛 equations for the power 𝑈 of each mode (where 𝑛 is the number of modes).

3. Experimental Validation of the Model

The experimental validation of the model is based on the measured characteristics of a 1.3 μm quantum dot (QD) oxide-confined (OC) VCSEL structure (Figure 1) as described in [24].

Figure 1 

A mesa structure of the 1.3 μm top-emitting oxide-confined (OC) In(Ga)As/GaAs quantum-dot GaAs-based VCSEL.

The current flow is defined by selective oxidation of an AlAs layer placed at antinode position of the mode. The active region is composed of five groups of three 8 nm thick In_0.15Ga_0.85As QWs, each of which containing one InAs QD layer. In each group, located close to the successive antinode positions of the optical standing wave within the cavity, the QWs are separated by 32 nm thick GaAs barriers. Additional single InGaAs QWs containing single QD layers are placed at the beginning and the end of the active region. The cavity is bounded by two Al_0.9Ga_0.1As/GaAs distributed Bragg reflectors (DBRs).

Figure 2 presents the comparison of the results obtained with the aid of the self-consistent model with measurement. Good agreement between theoretical and experimental results is evident. In the model we do not consider spontaneous emission, therefore the theoretical total emitted power increases sharply from the lasing threshold, which is slightly less than 6 mA. While the emitted power increases in a single-fundamental-mode regime, the carrier concentration decreases in the center of the active region, where stimulated recombination is the most intense because of spatial hole burning. Hence, the first-order mode LP₁₁ can be excited since more carriers are located closer to the radial borders of the active region.

Figure 2 

RT (18°C) CW light-current characteristic of the 10 μm diameter 1.3 μm InAs/GaAs quantum-dot GaAs-based oxide-confined (OC) VCSEL: crosses: experimental points, lines: theoretical curves; dash line: the LP01 mode, dot-dash line: the LP11 mode, and solid line: both the LP01 and the LP11 modes. Exactness of an output power measurement is not worse than 0.01 mW and that of an operation current is −0.1 mA.

It is clearly seen that the LP₁₁ mode appears at around 8 mA. For larger currents, although the total efficiency remains almost unchanged, the fundamental mode power decreases, while the higher-order mode quickly increases. Very exact agreement between experimental and theoretical plots shown in Figure 2 confirms validity of computer simulation.

Figure 3 presents a dependence of the RT CW lasing threshold currents of the 1.3 μm top-emitting In(Ga)As/GaAs QD GaAs-based VCSELs on diameters 𝑑_𝐴 of their active regions. If we take into consideration that the final shape of an oxide aperture is not a perfect circle, but it resembles rather an ellipse, which is followed by an exactness of the 𝑑_𝐴 determination not better than about 2 μm, experimental results agree quite well with theoretical ones. The lowest lasing threshold of 4.95 mA has been anticipated for 𝑑_𝐴 = 6 μm. In the case of smaller VCSELs, the radial size of the active region is too small to effectively confine the optical field. This causes a radical increase in its penetration into the passive areas around the central active region leading to a considerable increase in the lasing threshold observed in Figure 3. In VCSELs with larger active regions, current injection into the active region becomes increasingly nonuniform (current-crowding effect close to the active-region edge), which is followed by an excitation of higher-order transverse modes with higher thresholds than the fundamental mode.

Figure 3 

Dependence of the RT (18°C) CW threshold currents 𝐼th of the 1.3 μm top-emitting In(Ga)As/GaAs quantum-dot GaAs-based OC VCSELs on diameters 𝑑𝐴 of their active regions: crosses: experimental points and dash line: simulation curve. Exactness of the diameter 𝑑𝐴 seems to be not better than about 2 μm because a final shape of an oxide aperture is not a perfect circle, but it resembles rather an ellipse.

4. Comparison of Self-Consistent Models Assuming Scalar and Vectorial Optical Approaches

Simple scalar optical approaches based on the general assumption of the optical field within the laser cavity in the form of a plane wave [22] reveal very good accuracy in comparison with experiment as shown in the preceding section. However, modern designs of VCSELs are often equipped with relatively small active regions to reduce their lasing thresholds. The cavity size becomes comparable to the laser wavelength so that the above plane wave assumption is no longer valid. Then, the much more complicated vectorial optical approaches have to be used, which require solving six coupled equations for all components of both electric and magnetic fields. These algorithms are very involved, which increases the time taken to complete calculations by a factor of up to one hundred as compared to the scalar calculations. However, simplified scalar optical models sometimes give approximately correct results even beyond their confirmed range of validity. In such cases it may be unnecessary to use the full vectorial model.

We examine the limits of the scalar optical approach by comparing the threshold of a nitride VCSEL with more rigorous vectorial optical calculation. The reliability of the optical models used in the self-consistent analysis, that is, the plane wave admittance method and the effective frequency method, has been proved by comparison of their results obtained for the benchmark VCSEL structure from the COST-268 with those reported in [24]. The main conclusion drawn is that larger discrepancy occurs for the modal gain than for the emission wavelength and that the difference between them becomes larger for higher-order modes and smaller apertures where diffraction processes are essential [25].

Let us consider a possible structure of the GaN-based VCSELs (Figure 4) emitting at 400 nm. The structure and parameters used in the simulations are described in [26]. Our self-consistent simulation model is used to determine possible room-temperature CW performance characteristics of GaN-based VCSELs. We consider the fundamental LP₀₁ scalar mode and the fundamental vectorial HE₁₁ mode as well as the first-order LP₁₁ scalar mode and corresponding three vectorial HE₂₁, TE₀₁, and TM₀₁ first-order modes [27]. These three vectorial modes degenerate to the LP₁₁ mode in the scalar limit, but, as the difference between these three modes is rather small, we restricted our simulation to the HE₂₁ mode only.

Figure 4 

Structure of the modeled GaN-based VCSEL with the active region composing 5 quantum wells (MQW) and the tunnel junctions (TJ).

RT CW threshold current of a particular mode 𝐼_th (Figure 5(a)) increases superlinearly with increasing 𝑟_𝐴, mostly as a result of a steadily increasing active-region size. A possible impact of the active-region temperature (at threshold, its maximal value is increased from about 320 K for 𝑟_𝐴 = 2 μm to as high as 375 K for 𝑟_𝐴 = 8 μm, see Figure 5(b)) is ambiguous. It is followed by slightly lower and somewhat shifted material optical gain and some deterioration of other material properties but, at the same time, the radial thermal focusing leads to better optical confinement. Nevertheless, for too large active regions (i.e., for 𝑟_𝐴 > 8 μm), the VCSEL does not reach the lasing threshold at all, probably because of its excessively high active-region temperature.

(a)

(b)

(c)

Figure 5 

(a) RT CW threshold current of the particular modes 𝐼th determined with the aid of the fully self-consistent approach for two, the lowest-order scalar and vectorial transverse, modes versus the active-region radius 𝑟𝐴. (b) The highest active-region temperature 𝑇th in the threshold of the particular modes determined versus the radius of the active region 𝑟𝐴. (c) Radial profiles of the RT CW threshold carrier concentration 𝑛th(𝑟) determined for VCSELs equipped with active regions of various radii 𝑟𝐴.

For relatively large active regions (Figure 5(a)), the first-order vectorial HE₂₁ mode becomes the lowest-threshold one. This is a consequence of the changes in the carrier-density radial profile for larger active regions (see Figure 5(c)) with its maximum close to the active-region edge. In addition, the first-order vectorial mode exhibits lower lasing threshold than that of the scalar one. The impact of diffraction losses in these large-size devices is reduced.

For a VCSEL design equipped with a tunnel junction and an n-type radial-current-spreading layer (instead of the p-type one), both the above optical methods have been found to give accurate results for the fundamental transverse mode, which confirms the general opinion that the scalar effective frequency method yields reasonable results even beyond the range of its confirmed validity. Larger discrepancies are found for first-order transverse modes and for larger (or very small) active regions. For the same carrier concentration, the modal gain of the fundamental mode is usually somewhat higher (and the lasing threshold slightly lower), when it is determined using the vectorial approach, than that in the case of scalar one. The above is partly compensated by diffraction losses, not included in scalar approaches, which are more important for higher-order modes and smaller active regions. Therefore, contrary to the fundamental mode, modal gain determined for the same carrier concentration is higher (and the lasing threshold lower) for the first-order mode using the scalar optical approach than in the case of the vectorial one, especially in smaller devices. To summarize this comparison, simplified scalar optical approaches have been found to give approximately good results for fundamental transverse modes, but their accuracy is somewhat worse for first-order and higher-order modes as well as for large or very small VCSEL diameter.

5. Simulation of Oxide-Confined VCSELs

In this paragraph we consider various designs of oxide-confined (OC) VCSELs using our self-consistent model. We explore structures that produce single-mode operation for relatively large optical apertures, which enables high-power lasing.

A typical double-intracavity-contacted structure (with two ring contacts) of a 1.3 μm OC GaAs-based (GaIn)(NAs)/GaAs quantum-well VCSEL with a single oxide aperture is shown in Figure 6 [28]. A radial selective oxidation is proposed to create an oxide aperture of radius 𝑟_𝐴 which serves as both an electrical (to funnel current spreading from annular contact towards the central active region) and an optical aperture (to confine an optical field in a radial direction). In the standard VCSEL structure, both the active region and the oxide aperture are located at the antinode positions of the optical standing wave (Figure 6). So the VCSEL is working as an index-guided (IG) device with a strong radial built-in index confinement.

Figure 6 

A typical double-intracavity-contacted structure of standard GaAs-based oxide-confined (OC) VCSEL with a single oxide aperture located at an antinode position of the optical standing wave within the VCSEL 1.5 λ cavity.

Lasing thresholds of cavity modes depend on the overlap of their intensity profiles with the optical-gain distribution. The latter is directly associated with a radial profile of the threshold current density injected into the VCSEL active region (Figure 7), becoming in these double intracavity contacted VCSELs more and more nonuniform with increasing active-region radius 𝑟_𝐴 [25, 26]. Therefore, for larger 𝑟_𝐴 it is found that higher-order transverse modes become the lowest-threshold modes, which limits the desired Guassian mode operation to relatively small active region diameter. These conclusions are confirmed in Figure 8 which shows the wavelength and maximum threshold gain of the lowest-threshold transverse LP modes versus the radius of the active region. As expected, fundamental mode operation is possible in small VCSELs with an active-region radius not exceeding 3 μm. An increase in the active-region size creates an increasing order of the lowest-threshold transverse mode.

Figure 7 

Radial profiles of the threshold current densities 𝑗th of the indicated lowest-threshold LP modes determined for different radii 𝑟𝐴 of the active region of the index-guided double-intracavity-contacted 1.5 λ cavity GaInNAs/GaAs DQW VCSEL.

Figure 8 

Wavelengths λ of the indicated lowest-threshold cavity modes (dashed curve) and their maximal threshold optical gain 𝑔th,max (solid curve) in the index-guided 1.5 λ cavity double-intracavity-contacted GaInNAs/GaAs DQW VCSEL.

The radial built-in index guiding (IG) is the strongest for an oxide aperture localized at an antinode position of the standing optical wave within the VCSEL cavity. Let us analyze the impact of the reduced waveguiding, which happens with a gradual shift of the aperture towards a node position in a VCSEL with a large active region diameter of 20 μm (Figure 9). As one can see for relatively small shifts of the aperture, the radial built-in IG is still sufficient to confine the optical field with a very low threshold gain. However, when the aperture is shifted too far from the antinode position, the IG becomes ineffective and is replaced by the gain guiding (GG) mechanism, which requires much higher gain to reach the lasing action. As shown in Figure 9, a dramatic increase in the threshold gain from about 1750 cm⁻¹ to as high as 4900 cm⁻¹ occurs with a shift of the aperture. At the same time, instead of the high-order LP transverse modes, the fundamental LP₀₁ mode becomes the lowest-threshold mode. A similar behaviour may also be seen in Figure 10 for a VCSEL with a much smaller active-region diameter of only 4 μm. In this case for all aperture positions, the fundamental mode remains the lowest-threshold mode, although its threshold for the GG regime is nearly six times higher than that for the IG one.

Figure 9 

Maximal RT CW threshold optical gain 𝑔th,max (solid line) and wavelength λ (dashed line) of the indicated lowest-threshold LP modes for 20 μm diameter 1.5 λ cavity double-intracavity-contacted GaInNAs/GaAs DQW VCSEL as a function of the distance of the oxide aperture from its antinode position within a laser cavity.

Figure 10 

Maximal RT CW threshold optical gain 𝑔th,max (solid line) and wavelength λ (dashed line) of the indicated lowest-threshold LP01 mode determined for the 4 μm diameter 1.5 λ-cavity double intracavity contacted GaInNAs/GaAs DQW VCSEL as a function of the distance of the oxide aperture from its antinode position within a laser cavity.

As expected, the single-fundamental-mode operation is achieved in IG VCSELs only for relatively small active regions but their lasing threshold remains very low even for VCSELs with the largest active regions considered. By contrast in GG VCSELs, the single-fundamental-mode operation is ensured for all devices, but their lasing thresholds are considerably higher.

The standard OC VCSEL with two identical oxide apertures [27] is equipped with a 3-wavelength-long cavity with apertures positioned exactly at the standing-wave antinode positions on both sides on its active region. The layered structure is very similar to that of the 1.5 λ cavity VCSEL with a single aperture: both DBR mirrors and the active region are identical, and only the n-type and p-type spacers are readjusted to create the 3 λ cavity and to locate both apertures at antinode positions. The operation of such a VCSEL with two identical oxide apertures has been theoretically analyzed in [28]. As expected, the threshold gain turns out to be somewhat lower than that of the 1.5 λ cavity IG VCSEL, but the desired single-fundamental-mode operation is again confined to small active regions of diameters not exceeding 7 μm, which is consistent with experimental results reported by Degen et al. [29]. A further increase in the active-region size is followed by a gradual increase in the mode order of the lowest-threshold transverse mode.

Let us consider a new approach [30, 31] proposed to enhance the single mode operation in larger OC VCSELs with two oxide apertures: one of them is shifted to the node position, where it behaves as an electrical aperture only, whereas the second aperture remains at antinode position (the spacer thicknesses are properly adjusted). In this design the diameters of the two apertures may be changed independently giving an additional degree of freedom. In particular, the radial current spreading between the annular n-side and p-side contacts and the central active region and the radial optical waveguiding may be optimized independently.

Since the kinetics of the oxidation is strongly dependent on temperature but also on compositions and thicknesses of oxidized (AlGa)As layers [32], it is possible that apertures of different diameters may be produced in one oxidation process.

The proposed VCSEL structure is a radial equivalent of the separate-confinement-heterostructure (SCH) [33, 34], which is created in the direction perpendicular to layer boundaries (Figure 11). In the SCH active region, both the carrier confinement in the thin active region and the field confinement within the wider cavity are produced independently by two (inner and outer) heterojunction pairs. By analogy, the proposed VCSEL structure shown in Figure 12 is able to independently create an active-region current-density profile and a radial optical index guiding with the aid of two different oxide apertures is called the separate-confinement-oxidation (SCO) structure [32] (Figure 11).

Figure 11 

A comparison of the separate-confinement-heterostructure (SCH) and the separate-confinement-oxidation (SCO) structures.

Figure 12 

The considered structure of the 1.3 μm 3 λ cavity SCO GaInNAs/GaAs DQW VCSEL.

The impact of a steadily reduced radius 𝑟_𝐸 of the electrical aperture on the RT CW lasing threshold of a 20 μm diameter SCO VCSEL is shown in Figure 13. As expected, when the oxide apertures are of the same size (𝑟_𝐴=𝑟_𝐸), the high-order transverse LP₁₂ mode exhibits the lowest RT CW lasing threshold. A decrease of 𝑟_𝐸 creates a gradual decrease in the mode order, accompanied by a considerable increase of the lasing threshold (Figure 13) and an active-region temperature increase Δ𝑇_𝐴,max (Figure 14). For diameters of the electrical aperture smaller than 15 μm, the desired stable single-fundamental-mode operation is achieved. The threshold optical gain is initially reduced because of a better overlap between both the distributions of the mode intensity and the optical gain; for smaller electrical apertures, the threshold optical gain increases again due to increased losses in the passive areas. For decreasing 𝑟_𝐸 in Figure 13, the threshold current, 𝐼_th, reaches its maximum value at 𝑟_𝐸≈7.5 μm and then continuously reduces because of the reduction of the active-region area. The wavelength of the lowest-threshold transverse mode depends on its penetration into the low-refractive-index oxide layer, which is larger for the higher-order transverse modes than for the lowest-order LP₀₁ one. That is why the LP₀₁ fundamental mode exhibits the longest wavelength.

Figure 13 

Impact of the radius 𝑟𝐸 of the electrical aperture on the RT CW maximal threshold gain 𝑔th,max and threshold current 𝐼th of the indicated lowest-threshold transverse LP modes determined for the large-size 20 μm diameter SCO VCSEL.

Figure 14 

Impact of the radius 𝑟𝐸 of the electrical aperture on the RT CW maximal active-region temperature increase Δ𝑇𝐴,max and wavelength λ of the indicated lowest-threshold transverse modes of the large-size 20 μm diameter SCO VCSEL.

Our theoretical model of the proposed SCO VCSEL configuration is consistent with the experimental work by Samal et al. [35]. As in the SCO VCSEL, two apertures of different diameters on both sides of the VCSEL active region are utilized leading to dramatically improved mode selectivity. However, there are insufficient details (such as cavity length) given in this paper to compare this structure with ours. Assuming that the active region is located at the antinode position, the larger aperture is placed close to the next antinode position, so it behaves as both an electrical and optical aperture. However, the second aperture is located relatively far from the active region. Nevertheless, Samal et al. [35] have observed single-fundamental-mode operation of their VCSELs with two oxide apertures, whereas their VCSEL with only one oxide aperture located close to the antinode position emits the high-order LP₃₁ transverse mode, which is consistent with our model.

6. Simulation of Photonic Crystal VCSELs

In this section we present numerical analysis of PhC VCSELs, which is aimed primarily at designing a PhC pattern that assures single-mode operation and simultaneously reduction of the threshold current with respect to the VCSEL without the PhC. To fully reflect the complexity of the phenomena taking place in PhC VCSEL operation, a self-consistent model comprising three-dimensional, vectorial optical, thermal, electrical, and recombination submodels is used in the analysis.

For simulation purposes we have chosen an InP-based, 1300-nm AlInGaAs VCSEL with tunnel junction carrier injection described in [36]. The lattice constant of PhC and hole diameter equal 3.2 and 1.6 μm, respectively. Although the results presented here are obtained for this particular device, we believe that such phenomena as the dependences of threshold current on electrical aperture and etching depth are quite general and would show qualitatively the same trends for other PhC VCSELs. However, the optimized parameters cannot be used straightforwardly for PhC VCSELs based on different material systems mostly because of the different wavelength of emission.

In the analysis we consider the fundamental (HE₁₁ mode which relates to scalar mode LP₀₁) and the first-order mode (HE₂₁ mode being the vectorial analog to the LP₁₁ mode) and we determine their threshold characteristics. Specifically we calculate their threshold currents as a function of etching depth, which is presented in Figure 15. An increase of the etching depth in the range from 0 to 5 μm leads to an increase in both threshold currents, which reveals that shallow holes cause leakage of the light and therefore higher modal losses. As soon as the holes become deep enough to interact with the mode, the threshold currents 𝐼_th of both modes start to decrease with etching depth. Holes etched through the top DBR and into the cavity assure threshold current close to a minimum value. Further increase in the etch depth does not change significantly the threshold current. The existence of such a plateau is due to the maximum mode confinement for the etching through the cavity, where the intensity of the mode is the largest.

Figure 15 

Threshold currents of the HE11 and HE21 modes as a function of the etching depth. The borders of the cavity and successive ten pairs of top and bottom DBRs are assigned by vertical lines. The region of strong discrimination and low threshold is assigned by hatched rectangle.

The threshold current dependences on etching depth are qualitatively similar for the two modes: they display maxima at a certain etching depth. The HE₂₁ mode interacts weaker with the PhC than the HE₁₁ mode; because the HE₂₁ mode is less confined, more light of HE₂₁ is escaping to the air holes hence deeper holes are necessary to achieve similar confinement to that observed for the fundamental mode. This behavior causes the first-order mode to suffer large losses for an etching depth, which is sufficient to confine the fundamental mode and to lower its threshold current. By proper optimization of the PhC, such as etching depth and optical aperture diameter, it is possible to simultaneously achieve low threshold current and strong higher-order mode discrimination in PhC VCSELs [37, 38]. In the hatched region in Figure 15, the threshold current for the HE₁₁ mode is lower and simultaneously the differences of the threshold currents of the HE₂₁ and HE₁₁ modes are larger than in the case of VCSEL without PhC; this indicates the optimal range of etching depths assuring high modal discrimination.

One can determine the region of optimal lattice constant and etching depth of a PhC pattern as shown Figure 16. For the sake of comparison, these regions are calculated for a tunnel junction VCSEL and for a proton-implanted VCSEL. These VCSEL structures differ by the relatively uniform gain spatial distribution in the first case and the strong current crowding in the second one [38].

(a)

(b)

Figure 16 

Regions of low threshold and strong mode discrimination (gray field) mapped in the plane of etching depth and optical aperture 𝑅𝐴 for 𝑎/𝐿=0.5 in the case of VCSEL with tunnel junction (a) and proton implanted (b).

7. Simulation of Photonic Crystal VCSELs for Enhanced Polarization Control

For several years, PhCs have been successfully used for providing birefringence and dichroism in Photonic Crystal fibers, although there are not many works for such PhC application in VCSELs. So far, only one such successful experimental attempt is reported in [39], where photonic crystal with elliptical holes is applied to provide stable polarization lasing. In that work, the stabilization of polarization with over 20 dB polarization mode suppression ratio (PMSR) is achieved with hexagonal lattice PhC with air holes elongated along the K and M crystallographic directions. Quite importantly, the PhC effect is accompanied by an anisotropic current injection that according to our investigations is the main reason of the observed dichroism. If we neglect this effect, we can observe that for the structure presented in Figure 17 modal dichroism is negligible; however, a slight frequency separation appears, as shown in Figure 18. This separation can be used to promote one of the polarizations due to the better overlap with the gain profile.

Figure 17 

Schematic diagram of elliptical photonic crystal.

Figure 18 

Resonant frequency separation of the x- and y-polarized modes for the structure with elliptical holes as a function of the 𝑑𝑥/𝑑𝑦 ratio.

As can be seen from Figure 18, the hole structure with elliptical holes provides approximately 180 GHz frequency separation and no dichroism. From a practical point of view, this is not enough to ensure stable single polarization of the emitted light. Hence, we propose another design, based on dichroic optical fiber [40]. The hole structure presented in Figure 19 provides over 7% dichroism and more than 400 GHz of frequency splitting between the two orthogonally polarized fundamental modes.

Figure 19 

Schematic diagram of birefringent photonic crystal for light confinement in polarization-stable VCSEL.

The optical mode is confined by the photonic crystal pattern of Figure 19 with C6 symmetry broken by four holes with diameter two times larger than the other holes. The factor of two-hole diameter enhancement is important, as its decrease quickly reduces the birefringence and dichroism, while further increase is hard to achieve due to the resulting thin walls between holes. The computed birefringence and dichroism of this structure are presented in Figures 20 and 21. It turns out that a small hole in the middle plays crucial role for the polarization stabilization as it increases the overlap of the mode with the region where the large holes are located and thus the crystal symmetry is broken. Figure 22 shows electric field profiles for the two orthogonal polarizations with different sizes of this hole.

Figure 20 

Frequency separation between the two orthogonal linearly polarized fundamental VCSEL modes as a function of the middle hole diameter.

Figure 21 

Relative loss/threshold gain difference between the two orthogonal polarized VCSEL fundamental modes as a function of the middle hole diameter.

Figure 22 

Electrical field profiles of both polarizations for the birefringent PhC structure for two different sizes of the middle hole.

8. Conclusions

We have presented a detailed self-consistent way of modeling the various physical phenomena that occur within a VCSEL. The simulation comprises thermal, electrical, optical, and gain submodels. All of the submodels are combined in a self-consistent manner. Such an approach allows consideration of nonlinear effects occurring during the laser emission, such as spatial hole burning. We first showed the conformity of the numerical model with experimental results. In particular, we have shown that the model reproduces very precisely the experimental characteristics of both emitted power as a function of current and threshold current versus the diameter of oxide aperture.

Modern designs of VCSELs are often equipped with relatively small active regions to reduce their lasing thresholds, and, as result, diffraction effects become significant. In order to precisely analyze such devices, more accurate vectorial models are necessary since simple scalar models are applicable only for optical apertures several times larger than the wavelength of the cavity mode. Our comparative analysis showed that the scalar model is accurate for the fundamental mode and standard devices. Higher-order mode analysis, devices with narrow or very wide apertures—for which subtle thermal issues affecting waveguiding become prominent—should be analyzed with vectorial optical models. Our self-consistent model has been used to design oxide-confined VCSELs that exhibit stable, single-mode operation. The typical VCSEL structure with an oxide layer positioned at an antinode position provides strong optical confinement, which favors high-order modes for broader than 6 μm optical apertures.

The way to support single-mode operation for broader apertures is through weakening of the optical confinement. That can be realized by shifting the oxide layer away from the antinode position. Such approach assures single-mode operation, at the expense of a rapidly growing threshold current. The design comprises two oxidations: one layer in a field node position and another one in an antinode position. The first narrower oxide aperture funnels the current, while the second thicker oxide aperture confines the mode. Such a design provides stable single-mode operation up to electrical diameters as large as 15 μm.

Another analyzed method to achieve stable single-mode operation that additionally reduces the threshold current is the introduction of a photonic crystal pattern of holes into the VCSEL. We showed that precise numerical simulation can predict the photonic crystal parameters that lead to improving the VCSEL performance with respect to mode discrimination and high power simultaneously. To assure not only the single-mode but also single-polarization operation, additional patterning of the photonic crystal is necessary. We presented the modeling of a photonic crystal VCSEL with an irregular geometry of holes, which enables large polarization mode suppression ratio.

Numerical modeling is a very powerful tool for designing VCSELs and exploring new phenomena. It is much cheaper than the trial and error method and allows for predicting behaviors of devices that do not yet exist. However, the assumptions used in deriving the models, the VCSEL parameters, and boundary conditions should be very closely related to the experiment. The step that should be always done first for a VCSEL model is its calibration with the measurements. This might be a time-consuming process but it is necessary to keep the theoretical analysis attached to reality.

Acknowledgments

The authors wish to thank Professor K. D. Choquette for his kind assistance in editing the manuscript. This work is supported by Grant no. N N515 533338 from the Polish Ministry of Science and Higher Education (MNiSzW) and by the Swiss National Science Foundation (SNF) through Grant SCOPES IZ73ZO_128019, COST MP0702.