1. Introduction

Previous attempts to simulate vertical cavity surface emitting lasers (VCSELs) have often concentrated on the optical problem (see, for example, [1,2]), as the small cavities are expected to exhibit interesting optical properties. There also have been several sophisticated optical models coupled to rudimentary electronic models (for example [3]), which attempt to capture the essential physics of the operation of the laser with a minimum of overhead. Finally, there have been a few attempts to comprehensively model VCSELs (for example [4]) using precise optical and electronic models. However, previous electronic models have lacked the full dynamic capabilities of the edge emitting laser simulator MINILASE II, which are necessary for obtaining accurate modulation response curves [5]. Also, previous optical models have involved approximations which are not valid for all VCSELs [4]. We describe in this paper an ongoing effort to comprehensively model VCSELs through an extension of MINILASE. The extension will include a new optical mode solver, partially based on a similar solver previously used by Lin and Deppe [2]. This full vector three dimensional mode solver finds active rather than passive cavity modes, treats both index and gain guiding equally well, and doesn’t require that the cavity geometry be separable. A simple non-separable index guided structure is analyzed to demonstrate the capabilities of the mode solver.

2. Electronic simulation

In MINILASE, the solution of the two dimensional electronic equations that determine electron and hole supply to the active lasing region is accomplished along the lines of silicon device simulation [6] with the addition of specific provisions that describe transport over abrupt heterojunctions [7] and capture in the active quantum well region [5]. The quantum wells are coupled to the classical transport equations by Bethe’s thermionic emission theory and by the inclusion of ballistic regions [5]. Within the quantum wells we solve a Boltzmann transport equation in energy space and a rate equation for phonon heating [8]. The resulting electrical simulator has been described in detail in [9] and represents the core of the MINILASE II simulation package that has been extensively tested through simulations of edge emitting lasers and subsequent comparisons with experimental results [9]. It provides the gain function including all known nonlinear effects that arise from transport phenomena: spectral hole burning, carrier and phonon heating, and finite capture rates, for example. The electrical simulator is expected to be unchanged for VCSELs and provides us with the carrier densities and gain that are required by the optical solver that is described as the main topic of this paper.

3. Optical simulation

3.1 Theory

MINILASE solves the photon rate equation for the intensity of a given optical mode. The photon rate equation for a single mode is

where S is the mode intensity, G is the gain (due to stimulated emission) for the mode under consideration, τ_phot is the photon lifetime of the mode, and R_spon is the amount of spontaneous emission which couples into the mode. In order to find G, τ_phot , and R_spon for the mode under consideration, the optical modes must be obtained from the Maxwell equations without field-independent sources (spontaneous emission). In the conventional semiconductor laser theory [10], a set of simplifying assumptions are made in order to find the optical modes. The cavity geometry is usually assumed to be separable to some degree, and the gain is often treated as a small perturbation to the permittivity ∊(r⃗, ω). In practice, the optical cavity is often closed in order to find the modes, and the effect of the gain and cavity loss on the spatial mode pattern is neglected altogether [9]. The mode patterns calculated in this way are known as the passive or cold cavity modes. The assumption that the gain distribution has little effect on the lasing mode pattern is well justified for strongly index confined devices [10]. However, many VCSEL devices may be strongly gain guided. Therefore it is desirable to include the gain in ∊(r⃗, ω) and find instead the active cavity modes of the VCSEL structure.

In order to solve Maxwell’s equations for the active cavity modes, we first separate ∊(r⃗,ω) into two parts as

where ∊_cav(r⃗, ω) is the potentially complex relative permittivity representing the passive cavity, and χ_gain (r⃗, ω) is the necessarily complex susceptibility representing the gain region of the laser. A susceptibility rather than a permittivity is used to represent the gain in order to clearly separate it from ∊_cav , for reasons that will become apparent shortly. We can write the vector Helmholtz equation for the electric field as

where J⃗_stim,(r⃗,ω) = -iω ∊_o χ_gain (r⃗,ω)E⃗(r⃗,ω) is the field-dependent current representing stimulated emission in the gain region. The solution of the inhomogeneous wave equation Eq. (3) can be represented as

where G⃡_cav (r⃗,r⃗′,ω) is the tensor Green’s function for radiation from a current source within the open VCSEL cavity defined by ∊_cav .

Eq. (4) can be viewed as a condition that the gain susceptibility χ_gain must satisfy in order that the corresponding mode E⃗ lase (be self-supporting) in the absence of spontaneous emission. With spontaneous emission present, the gain (and the mode) must approach this condition asymptotically. In order to arrive at a more useful form of Eq. (4), we note that spatial distribution of the gain susceptibility χ_gain can be calculated from the electronic system using MINILASE. We therefore write χ^gain (r⃗,ω) = K(ω)${\chi }_{\mathit{\text{gain}}}^{\left(0\right)}$(r⃗,ω), where ${\chi }_{\mathit{\text{gain}}}^{\left(0\right)}$ is the normalized gain susceptibility pattern, and K is a complex constant that allows χ_gain to meet the lasing condition given as Eq. (4). Now Eq. (4) becomes

where G⃡ (ω) is an integral operator defined by

Eq. (5) is now an eigenvalue problem; the eigenfunctions E⃡ are the mode field patterns and the eigenvalues K are the gain amplitudes which define the lasing condition for the corresponding modes. This eigenmode equation is unusual as the eigenvalues are complex gain eigenvalues rather than complex frequency eigenvalues.

The complex gain eigenvalues K are a continuous function of frequency. In order to obtain the frequencies at which the VCSEL can actually lase, the complex gain susceptibility necessary for the optical system must be matched to the complex gain susceptibility available from the laser electronic system in both magnitude and phase. This match will only occur at discrete frequencies. This determines the lasing frequencies of the laser. In addition, the required photon rate equation parameters may be calculated for each mode obtained in this fashion through use of the lasing mode pattern E⃗ and the gain susceptibility required to lase χ_gain . For example, we may obtain the photon lifetime by using the fact that gain equals loss at threshold; hence 1/τ_phot will equal the total gain of the active cavity mode at threshold, which may be calculated as the overlap integral of the mode’s normalized field intensity pattern times the gain susceptibility required to lase.

It is possible to repartition the VCSEL susceptiblility to take part of ∊_cav out of the Green’s function and put it into the eigenvalue problem instead. Suppose we have a localized passive cavity feature which we don’t want to include in the passive cavity Green’s function G⃡_cav We may then split the electric permittivity further as

where the new term χ_cav represents a localized deviation from the cavity structure defined by ∊_cav . This susceptibility then turns up in Eq. (3) as a second source term, and Eq. (5) now takes the form of a generalized eigenvalue equation:

There are two primary reasons to use this generalized formulation. First, if some localized regions of the laser cavity change with increasing bias (e.g., localized temperature related index changes) it would be useful to put these changes into the eigenvalue problem rather than the Green’s function, to avoid the need to recalculate Green’s functions. Second, there are a wide range of VCSELs which can be effectively modeled as very simple cavity structures containing localized deviations from the simple structure. The simple cavity structure will have an easily calculable Green’s function. The localized deviations from the simple structure can be partitioned into χ_cav Thus we greatly simplify the task of finding the Green’s function with only a modest corresponding increase in the difficulty of the eigenvalue problem. Making use of this allows greater computational speed and may in fact produce greater accuracy as well, by avoiding the complicated and slow numerical techniques needed to find the tensor Green’s function for a generalized cavity structure.

 

Fig. 1. Index guided VCSEL structure used to generate results described below. λ_o =1 μm

Download Full Size | PDF

3.2 Implementation

We present here some sample results from our current implementation of the gain eigenvalue method. Although these simulations did not include the electronic solver, useful information and insight can be gained by examining the optical modes and corresponding gain eigenvalues.

We used the gain eigenvalue method to find the optical modes of a simple but non-separable problem, the index guided VCSEL diagrammed in Fig. 1. The underlying cavity structure ∊_cav is planar, and therefore an analytic Green’s function is available [11]. The wavelength-tall optical cavity is bounded above and below by distributed Bragg reflectors (DBRs). There is an oxide layer surrounding the active region in the plane of the active region (similar to the model of Lin and Deppe [2]). The oxide layer is slightly lossy in order to enforce the radiation condition as well as to physically model optical losses in the VCSEL structure. For our initial simulations, we set the normalized gain susceptibility pattern ${\chi }_{\mathit{\text{gain}}}^{\left(0\right)}$ to be uniform in a rectangular active region, and zero elsewhere. Thus for later convenience we set ${\chi }_{\mathit{\text{gain}}}^{\left(0\right)}$(r⃗) ≡ -i within the 9μm ×10μm×3nm quantum well active region, and ${\chi }_{\mathit{\text{gain}}}^{\left(0\right)}$(r⃗) ≡ 0 elsewhere. The entire layer which contains the oxide and the active region is given the background index of the oxide, which is effectively raised to the index of GaAs inside the active region by setting χ_cav inside the active region to be equal to the difference between the two permittivities.

Once the Green’s function for sources within the active region was calculated, the gain eigenvalue problem of Eq. (8) was discretized and solved. Using a moment method, the field was represented as a sum of known expansion functions with unknown coefficients, and the inner product of both sides of Eq. (8) was taken with a set of known testing functions. The resulting dense matrix eigenvalue problem was solved using a commercially available eigenvalue solver (IMSL).

Figure 2 shows the real and imaginary parts of the first six gain eigenvalues K(ω) for this structure (the first two are essentially degenerate) plotted at constant intervals in frequency. Approximating the gain susceptibility available from the electronic system of the laser under bias as purely imaginary (neglecting the Kramers-Kronig relation for now), lasing is only possible at frequencies at which one of the gain eigenvalue curves crosses the real axis, such that χ_gain (r⃗,ω) = K(ω)${\chi }_{\mathit{\text{gain}}}^{\left(0\right)}$(r⃗) is also purely imaginary. The gain eigenvalue at the crossing frequency defines the gain susceptibility required to lase. The large area active region considered in this example produces a very tight spacing of lasing modes in frequency, but the gain amplitudes required for lasing are significantly different because of the different overlaps of the self-consistently calculated lasing modes with the active region. Figure 3 shows the field patterns of the first and third lasing modes inside the active region. These modes are full-vector modes of the three dimensional cavity.

 

Fig. 2. First six complex gain eigenvalues plotted at discrete intervals in frequecy; the first two are essentially degenerate. The gain eigenvalues corresponding to the first three modes are circled and labeled. A mode can lase only at the frequency at which its gain eigenvalue curve crosses the real axis. Therefore modes 1&2 lase at λ ≃ 1.0007λ_o and mode 3 lases at λ ≃ 0.9999λ_o . The magnitude of K at the crossing frequency determines the amount of gain required for the corresponding mode to lase. Thus we see that mode 3 requires more gain to lase than modes 1&2.

Download Full Size | PDF

 

Fig. 3. Transverse mode patterns of the (a) first and second, and (b) third optical modes. The first mode is primarily x̂-polarized and the second is primarily ŷ polarized. The third mode is primarily ŷ-polarized. The next three modes shown in Fig. 2 are permutations of the third mode, with either field pattern or polarization rotated 90 degrees, or both.

Download Full Size | PDF