## Abstract

An approximation to the Maxwell-Semiconductor Bloch equations is used to model transverse mode dynamics of vertical-cavity surface-emitting lasers (VCSELs). The time-evolution of the spatial profiles of the laser field and carrier density is solved by a finite-difference algorithm. The algorithm is fairly general; it can handle devices of any shape, which are either gain or index guided or both. Also there is no a priori assumption about the type or number of modes. The physical modeling includes the nonlinear carrier dependence of the optical gain and refractive index and dispersion effects on the gain and the refractive index are also included.

© Optical Society of America

## 1. Introduction

VCSELs have become increasingly important for many applications, such as optical interconnects and optical data storages. In many applications, high beam quality or high output powers or both are critical requirements. To help the design and optimization of VCSELs for such applications, it is essential to have a simulation model that includes the time evolution of the transverse space dependence of the beam. Such models allow arbitrary transverse profiles of laser intensity or carrier density to develop in time without any a priori assumption about the number and types of transverse modes. In addition, extensive research [1] in the past has shown that many-body effects are important in influencing the optical gain and the refractive index of semiconductor quantum wells. It is thus important to include the key many-body effects in simulation models for semiconductor lasers.

Current methods of VCSEL simulation are limited in various ways from this more complete modeling. For example, some methods select a few transverse modes beforehand and then solve their time evolution either by ordinary or partial differential equations [2]. Other methods solve an eigenvalue problem for the Helmholtz equation, which is uncoupled to the material equations. Finally, other methods solve time independent, coupled rate equation models, which contain diffractive terms (in the wave equation) and diffusive transverse terms (in the carrier density equation). Common limitations to most of these approaches are that the gain and refractive index are linearized with respect to the density (with the use of phenomenological coefficients), the dispersions with respect to wavelength of the gain and refractive index are neglected and finally that the many-body effects are not included.

There are two main methods to developing simulations with the complete modeling of the time evolution of arbitrary transverse profiles, as mentioned above. The first method is to solve the coupled Maxwell-Semiconductor Bloch equations [1] directly. This approach has been employed either on an envelope equation approximation or on the original Maxwell-Semiconductor Bloch equations [3–5]. Due to the time scales and the large set of equations involved, it is only feasible to solve this set of equations for a time development of up to picoseconds for a 1-dimensional simulation. Current computer capabilities render it impractical to simulate the time evolution of the transverse profiles of semiconductor lasers for a time development on a scale of several nanoseconds, and it is on this time scale that most of the interesting features of lasers occur.

The second method solves an approximation to the Maxwell-Semiconductor Bloch equations. These approximate equations are called the Maxwell Effective Bloch Equations [6,7] This model combines the space-time resolution with microscopic many-body effects in a way that is both accurate enough and computationally manageable, and it has been used to simulate edge-emitting high power lasers [8,9]. In this approach, partial differential equations are solved in space and time for the coupled wave and material functions. Both the nonlinear carrier density dependence and the dispersions with respect to wavelength of the gain and refractive index are included in the modeling. This bottom-up approach uses measured material parameters and quantum well structure parameters, with the number of free parameters minimized. This simulation can resolve the time evolution down to picoseconds of the transverse profiles. Finally, this simulation is versatile in being able to readily model devices with any shape of active region, which are either gain or index guided, and without any assumption of the type or number of transverse modes.

In this paper, this model is applied to VCSELs. The Maxwell-Effective Bloch Equations model is adapted to VCSEL simulations by assuming a fixed mode profile along the longitudinal direction in VCSELs. This approximation is valid for VCSELs because after about a picosecond, the longitudinal mode structure is established [10,11]. There is not very much deviation from the assumed mode structure thereafter, because of the very large mode spacing between the adjacent longitudinal modes. In these simulations, the VCSELs are based on InGaAs/GaAs quantum well structures.

A finite-difference algorithm has been developed for solving the Maxwell Effective Bloch Equations [12,13]. Preliminary calculations [12] have been performed for two current apertures; a circular aperture with diameter 7.5 microns and a square aperture with sides of length 7.5 microns also. For each aperture three levels of pumping current were applied to obtain steady light fields of the fundamental modes and the light-current characteristics near threshold were determined. Next calculations [13] were performed for circular apertures at 5.0, 7.5, and 10.0 microns in diameter and at various pumping levels. The resulting light fields were time dependent with characteristic time scales on the picosecond time scale. Different types of transverse mode dynamics occurred that included multiple mode beatings, rotating waves, convective waves, and nonlinear locking of transverse modes.

In this paper [14], the effects of including index confinement and using ring shaped current apertures are studied. In these calculations, it was found that the average light output on the nanosecond time scale was steady, even though on the picosecond time scale the light output was very time dependent and in addition at relatively high current pumping levels, the light output appeared chaotic-like. Movie clips [15] show the transverse mode dynamics of the cases studied in this paper.

References 16 and 17 discuss experimental results concerning transverse mode dynamics of VCSELs.

## 2. Governing Equations

After assuming a fixed mode profile along the longitudinal direction in VCSELs, the resulting Maxwell Effective Bloch Equations model is given by the following equations [6,7,10,11].

Here *E* is the complex laser field envelope amplitude,
*N* is the total carrier density *P*
_{0} and
*P*
_{1} are the effective material polarization functions
that have been constructed from microscopic theory [6], *δn*(*x, y*) is
the guiding index profile, ${\nabla}_{\perp}^{2}$ is the
Laplacian in the transverse plane, (x,y), c is the speed of light,
*∊*
_{0} is the permittivity of free space, i
is the complex number √-1, *ω*_{c}
is
the optical carrier wave frequency in radians per seconds,
*n*_{b}
is the background index of refraction,
*n*_{g}
is the group index of refraction,
*∊*_{b}
=${n}_{b}^{2}$ is the background relative permittivity,
K=*ω*_{c}*n*_{b}*/c* is the optical
wavenumber in the cavity with a background index of refraction
*n*_{b}*, k* is the cavity loss,
*D*_{N}
is the carrier diffusion coefficient, γ_{n} is the nonradiative decay constant or carrier loss rate due to spontaneous
and nonradiative processes, *η* is the quantum efficiency,
e is the electron charge,
*ħ*=*h*/2*π*,
where h is Plank’s constant, Γ is the confinement factor, and
*L* is the cavity length.

Many-body effects and quantum well structure information are contained in the density
dependent coefficients _{X0}(*N*), the
effective background susceptibility, with real and imaginary parts
_{X0,r}(*N*) and
_{X0,i}(*N*)
respectively; Γ_{1}(*N*), the gain bandwidth,
*ω*
_{1}(*N*), the detuning, and
*A*
_{1}(*N*), the strength of the
Lorentzian oscillator. The theoretical basis for the Maxwell-Effective Bloch
Equations and their derivation from the semiconductor Bloch equations and
Maxwell’s equations is given in Reference 6. Also, the derivation of the
five density dependent coefficients, which model the optical sus-ceptibility,
_{X}(*ω,N*), is given in Reference 6.

#### 2.1 Optical Susceptibility

The five density dependent coefficients that describe the optical susceptibility in the Maxwell-Effective Bloch Equations are determined in two steps. In the first step, the optical gain and refractive index are computed using the specifications of the material compositions and other material parameters, including the microscopic many-body interactions, the bandstructure for the valence bands, and the specifics of the quantum well heterostructure.

Using current standard procedures, [1] for a given total carrier density, N, the semiconductor
Bloch equations are solved for the distribution functions,
${n}_{k}^{e}$
(*t*) and
${n}_{k}^{h}$
(*t*), the
occupation probabilities for electron and holes respectively, and
*p*_{k}
(*t*), the interband dipole
expectation function with momentum *ħk*. Next,
*P*(*t*), the induced polarization is computed [1] from the
*p*_{k}
(*t*). Then the optical
susceptibility _{X}(*ω;N*) is determined
from the Fourier transforms of the electric field,
*E*(*ω*), and the induced
polarization, *P*(*ω*) by the equation

From the optical susceptibility, the refractive index,
*δn*(*ω;N*), and the
optical gain, *G*(*ω;N*), are
determined by the equation,

In the second step, the optical susceptibility function,
_{X}(*ω,N*) is approximated by a
background susceptibility _{X0}(*N*) and one
Lorentzian oscillator _{X1}(*ω;N*).

Now the five density dependent coefficients have been determined for one value of the total carrier density N. This two step procedure is repeated for several values, (approximately twenty five values in the present case), of N to determine these functions over the range of interest. In this derivation, plasma heating is ignored, the approximations are made to locally fit the gain peak for each value of N and the distribution functions are assumed to be quasi-Fermi in solving the semiconductor Bloch equations. By using this method, the five material functions were determined. The resulting gain curves obtained by these Lorentzian approximations are shown in Fig. 1 for four values of the total carrier density N, (depicted by solid lines), and they are compared to the gain curves obtained by solving the semiconductor Bloch equations, (depicted by dashed lines), for which the material is described by an InGaAs/GaAs quantum well structure.

## 3. Computed Results

For these calculations, some of the parameters in the governing equations are the
following. The two parameters that seem to be the most important in determining the
dynamics were Γ_{1}(*N*), the gain bandwidth,
1/Γ_{1}(*N*) was approximately 13.8
femtoseconds initially in the current aperture area and the light dynamics changed
on that time scale and γ_{n}, the nonradiative decay constant, 1/γ_{n} was 1 nanosecond, and at that time scale the laser field was approximately
steady. The light wave length was 0.98 microns, the circular current apertures was
10.0 microns in diameter, the cavity length *L* was 144 nanometers
and the confinement factor was one fourth. The value used for index guiding was
-0.05, which was the decrease in the refractive index outside of the active region.

For each case, a movie clip [15] can be viewed to see the transverse mode dynamics.

#### 3.1 Case 1. Disk without index guiding

The first case was of a disk shaped current aperture without index guiding and the pumping current was 0.55 milliamps. Figure 2 shows the resulting steady laser field. The field rose to a fairly sharp high peak. The most intense light is the lightest in the figure and the lower intensity is darker. Outside the darker area, there is no light.

#### 3.2 Case 2. Disk with index guiding and low pumping current

The second case was of a disk shaped current aperture with index guiding and the pumping current was again 0.55 milliamps. Figure 3 shows the resulting laser field at one instant of time. The index guiding has a confining effect on the light field and the more intense field becomes dynamic with the peak making a generally rotating motion and the peak amplitude heaving to larger and smaller values. Figure 4 shows the average light field over 3.3 nanoseconds. Now the bright center is broader than in Case 1, but the peak value is lower.

#### 3.3 Case 3. Ring with index guiding and low pumping current

The third case was of a ring shaped current aperture with index guiding. The outer radius was again 10 microns, the inner radius was 4 microns and the pumping current was 0.48 milliamps. Figure 5 shows the resulting laser field at one instant of time. There are now 4 peaks of intensity located symmetrically around the ring current pumping area. Again, the field is dynamic with the 4 peaks generally rotating around the ring and the 4 peak amplitudes heaving to larger and smaller values in unison. Figure 6 shows the average light field over 3.3 nanoseconds. Now there is a bright ring shaped area over the current aperture with a dark center.

#### 3.4 Case 4. Disk with index guiding and high pumping current

The fourth case was of a disk shaped current aperture with index guiding and the pumping current was now raised to 2.35 milliamps. Figure 7 shows the resulting laser field at one instant of time. Now there are multiple peaks moving about in a chaotic-like manner, (but with a general rotation of the peaks on the rim), with some peaks disappearing and others appearing and their amplitudes rising and falling. Figure 8 shows the average light field over 3.3 nanoseconds. In comparison to Case 2, which had a lower pumping current, here the bright center of the average field has an even brighter circular edge, outside of which the light field dies away quickly.

#### 3.5 Case 5. Ring with index guiding and high pumping current

The final case was the same ring shaped current aperture of case 3 but the pumping current was now raised to 2 milliamps. Figure 9 shows the resulting laser field at one instant of time. There are now many peaks of intensity located around the ring current pumping area. They have various amplitudes and are generally slightly rotating, (i.e. slightly rocking back and forth), around the ring with some peaks disappearing and others appearing and their amplitudes rising and falling. Figure 10 shows the average light field over 3.3 nanoseconds. As in Case 3, there is a bright ring shaped area over the current aperture with a dark center.

## 4. Conclusion

A finite-difference algorithm has been developed for solving the Maxwell Effective Bloch Equations. This algorithm was applied to VCSELs. Calculations were performed for circular current apertures that were 10.0 microns in diameter without and with index guiding and both for disk and ring shaped current injection areas. Various steady current pumping levels were used. For many cases, the resulting light fields were time dependent with characteristic time scales on the picosecond time scale. Different types of transverse mode dynamics occurred that could be described in many cases as multiple moving and heaving peaks, although in the case of the ring shaped current aperture, the multiple peaks of light intensity tended to rotate around the ring while heaving in amplitude. At relatively high current pumping levels, the dynamics of the light output appeared chaotic-like. However, in all cases, the average light output on the nanosecond time scale was essentially steady.

## References and links

**1. **W. W. Chow, S. W. Koch, and M. Sargent, *Semiconductor Laser Physics*, (Springer, Heidelberg, Berlin, 1994). [CrossRef]

**2. **J. Y. Law, G. H. M. van Tartwijk, and G. P. Agrawal, “Effects of transverse-mode competition on the injection dynamics of vertical-cavity surface-emitting lasers,” Quantum Semiclass. Opt. , **9**, 737 (1997). [CrossRef]

**3. **P. M. Goorjian and G. P. Agrawal, “Computational Modeling of Ultrashort Optical Pulse Propagation in Nonlinear Optical Materials,” Paper NME31, *Nonlinear Optics: Materials, Fundamentals and Applications*, 11, 1996 OSA Technical Digest Series, Washington, D.C., 1996, 132–133.

**4. **P. M. Goorjian and G. P. Agrawal, “Computational Modeling of Ultrafast Optical Pulse Propagation in Semiconductor Materials,” Paper QThE9, Quantum Optoelectronics, Spring Topical Meeting, OSA, Washington, D. C, Nevada, March 17–21, 1997.

**5. **P. M. Goorjian and G. P. Agrawal, “Maxwell-Bloch Equations Modeling of Ultrashort Optical Pulse Propagation in Semiconductor Materials,” Paper WB2, OSA 1997 Annual Meeting, Washington, D. C, October 12–17, 1997.

**6. **C. Z. Ning, R. A. Indik, and J. V. Moloney, “Effective Bloch-equations for semiconductor lasers and ampliers,” IEEE J. Quantum Electron. **33**, 1543 (1997). [CrossRef]

**7. **T. Rossler, R. A. Indik, G. K. Harkness, J. V. Moloney, and C. Z. Ning, “Modeling the interplay of thermal effects and transverse mode behavior in native-oxide-confined vertical-cavity surface-emitting lasers,” Phys. Rev. A **58**, 3279 (1998). [CrossRef]

**8. **C. Z. Ning, J. V. Moloney, and R. A. Indik, “A first-principles fully space-time resolved model of a semiconductor laser,” Quantum Semiclass. Opt. , **9**, 681(1997). [CrossRef]

**9. **A. Egan, C. Z. Ning, J. V. Moloney, R. A. Indik, M. W. Wright, D. J. Bossert, and J. G. McInerney, “Dynamic Instabilities in MFA-MOPA Semiconductor Lasers,” IEEE J. Quantum Electron. **34**, 166, (1998). [CrossRef]

**10. **C. Z. Ning, S. Bischoff, S. W. Koch, G. K. Harkness, J. V. Moloney, and W. W. Chow “Micro-scopic Modeling of VCSELs: Many-body interaction, plasma heating, and transverse dynamics,” Optical Engineering, April, 1998.

**11. **C. Z. Ning, R. A. Indik, and J. V. Moloney, “A self-consistent approach to thermal e ects in vertical-cavity surface-emitting lasers,” J. Opt. Soc. Am. B **12**, 1993–2004, 1995. [CrossRef]

**12. **P. M. Goorjian and C. Z. Ning, “Computational Modeling of Vertical-Cavity Surface-Emitting Lasers,” Paper Thc15, Nonlinear Optics Topical Meeting, Kauai, HI, August 9–14, 1998.

**13. **P. M. Goorjian and C. Z. Ning, “Simulation of Transverse Modes in Vertical-Cavity Surface-Emitting Lasers,” 1998 Annual Meeting of the Optical Society of America, Washington, D. C, October 5–9, 1998.

**14. **P. M. Goorjian and C. Z. Ning, “Transverse Mode Dynamics of VCSELs through Space-Time Simulation,” Paper 3625–45, Integrated Optoelectronic Devices, Photonics West, 1999, (SPIE), San Jose, CA, January 23–29, 1999.

**15. **P. M. Goorjian and C. Z. Ning, “Transverse Mode Dynamics of VCSELs through Space-Time Simulation,” http://science.nas.nasa.gov/egoorjian/Pub/pub.html

**16. **C. J. Chang-Hasnain, J. P. Harbison, G. Hasnain, A. C. Von Lehmen, L. T. Florez, and N. G. Stoffel, “Dynamic, polarization, and transverse mode characteristics of vertical cavity surface emitting lasers,” IEEE J. Quantum Electron. **27**, 1402–1409, (1991). [CrossRef]

**17. **Y. Satuby and M. Orenstein, “Small-Signal Modulation of MultitransverseModes Vertical-Cavity Surface-Emitting Lasers,” IEEE Photonics Tech. Letters , **10**, 757–759, (1998). [CrossRef]