In this example we demonstrate a workflow for simulating the L-I curve of an InGaAsP-InP multiple-quantum well (MQW) edge emitting laser (EEL). Calculated LI curves and lasing frequencies at different temperatures are compared to results reported in literature.

## Overview

Understand the simulation workflow and key results

Laser library and MQW module are used to model a Fabry-Perot InGaAsP-InP MQW ridge laser using scripts. L-I curve and spectrum of the laser are calculated at different temperatures. The 1D TWLM laser model is suitable for simulating lasers with traveling wave geometry, such as edge emitting lasers. The MQW solver is suitable for simulating material gain in III-V semiconductors with zincblende crystal structure. The simulations are isothermal, suitable for lower power or pulsed lasers without pronounced self-heating at high injection currents. Simulations can be performed at different temperatures. Links to laser examples with other types of cavities, such as DFB and DBR, are included as well as a list of other laser result types not shown here.

### Step 1: Optical Mode Simulation

The first step is to calculate the optical mode profile and extract the effective index and group index of the fundamental (TE) mode as well as its confinement factor with respect to the gain medium. These are all calculated in the vicinity of the nominal target lasing frequency and are expected to be locally weakly frequency dependent. This calculation can be done using the FDE solver in MODE.

### Step 2: Gain Medium Simulation

Using the MQW Gain Solver, a 4x4 k•p calculation of the electronic bandstructure in the MQW gain medium is performed and the electronic bandgap, stimulated, and spontaneous emission spectra are extracted as a function of carrier density and temperature. MQW can be run from any product. In this case it can be run either in FDE (needed in the previous step) or INTERCONNECT (needed for the following step). Due to high barrier thickness, QWs can be considered uncoupled. In this case calculations can be performed for a single QW to reduce simulation time. At the end of calculations, the results are scaled to represent the full number of QWs. The electric field is set to zero, which is a good approximation for the forward bias laser diode near and above threshold.

### Step 3: 1D Traveling-Wave Laser Simulation

Using the TWLM element running in INTERCONNECT a time-domain 1D laser simulation is performed using a sweep over different drive currents. The optical power emitted in steady state may then plotted as a function of drive current to generate the L-I curve. The spectrum can be plotted at different times.

## Run and results

Instructions for running the model and discussion of key results

### Step 1: Optical mode calculation

- Launch MODE.
- Open and run the script file
*input.lsf*. This will initialize input data and simulation configuration parameters that are required for all the three steps. This data will be saved in*Piprek2000OQE.json*which will be called later on by scripts in different steps when needed. - Open and run the script file
*runFDE.lsf.*This will calculate the effective index as a function of frequency and then the group index using a finite difference approximation.

The effective and group indices, as well as the confinement factor used in the subsequent steps are stored in * wgSweep.neff, * * wgSweep.ng * , and * wgSweep.confinement_factor * , respectively. These results will also be stored in the file * Piprek2000OQE_opt.json * .

### Step 2. Gain medium calculation

- While still in MODE, open and run the script file runMQW.lsf to calculate electronic bandgap, stimulated, and spontaneous emission spectra as a function of carrier density and temperature. Alternatively, this script can be run from INTERCONNECT.

This script first loads the results from the eigenmode simulation if necessary, as the effective index will be required for the ensuing calculations. It will then set up the layer structure, calculate the material properties for the MQW gain region using the buildmqwmaterial script command and customize some of the properties using the values from the reference, and then run the MQW solver's mqwgain script command to calculate the gain curves, spontaneous emission spectra, and bandgaps as a function of temperature and the average carrier density in the gain region. These results are also stored in a Lumerical json file * Piprek2000OQE_mqw.json * .

Following plots represent spontaneous emission and stimulated emission vs frequency at different carrier densities are generated by the script (at 293K).

### Step 3. 1D Traveling-Wave Laser Simulation

- Launch INTERCONNECT.
- Open and run the script file
*runTWLM.lsf*. - Right clicking on the sweep, selecting power, and then visualize. This is the L-I curve for the laser operating at 293K (20C).
- Right clicking on the sweep, selecting spectrum, and then visualize. For a desired current, spectrum of the laser can be plotted at 293K (20C).
- Open the
*input.lsf*and set the parameter*twlm.temperature*to 333 to produce the L-I curve at 333K (60C). - Find the parameter twlm.li_currents and uncomment the list of currents more appropriate for the producing the L-I curve at 333K.
- Run the input.lsf file again.
- Run the runTWLM.lsf file.
- Right clicking on the sweep, selecting power, and then visualize. This is the L-I curve for the laser operating at 333K (60C).
- Right clicking on the sweep, selecting spectrum, and then visualize. For a desired current, spectrum of the laser can be plotted at 333K (60C).

The runTWLM.lsf script loads the results from the previous two simulations as well as the INTERCONNECT template file. This script will calculate the radiative and non-radiative recombination rates, create recombination files and load those parameters into the TWLM. It will also populate the other required parameters for the TWLM element, electrical GAIN element (used to model current leakage), root element, and parameter sweep. It will then run a very short simulation in order to allow the TWLM element to perform the gain and spontaneous emission fitting. Once the fits are performed the results are stored in the TWLM to prevent refitting them in each iteration of the sweep. The diagnostic fit results are also saved as datasets in the file called * Piprek2000OQE_293.ldf * . Finally, the sweep is launched.

Since the sweep can be parallelized, in order to shorten the simulation time it is possible to configure resources in Simulation menu to use multiple cores configuration. The easiest way is to simply add several resources (e.g. 3-4), depending how many cores are available, with default options (localhost for your local machine and 1 process per resource). Using more than 1 process per resource (multi-threading) may not always produce the desired speed-up.

Figure below shows the laser spectrum in steady state at 293K with a driving current of 0.2 A.

Figure below shows a comparison between simulated L-I curve at two temperatures and results reported in Reference 1. Also the threshold currents, slopes and peak wavelengths from the present Lumerical simulation as well as those from Reference 1 are presented in table below.

Lumerical | Reference 1 | |

Threshold Current [A] at 293 K | 0.12 | 0.114 |

Slop [W/A] at 293K | 0.26 | 0.24 |

Peak Wavelength [nm] at 293K | 1527 | 1536 |

Threshold Current [A] at 333 K | 0.245 | 0.247 |

Slop [W/A] at 333K | 0.24 | 0.19 |

Peak Wavelength [nm] at 333K | 1556 | 1552 |

## Important model settings

Description of important objects and settings used in this model

Below is a subset of important model settings. Information about other settings can be found in the comments in input.lsf script file.

### Multi-quantum-well region in FDE

Because the overall FDE domain is much larger than the MQW thickness we include the MQW domain as a single chunk of material with refractive index averaged over different quantum wells and barriers. This makes the simulation much faster at a reasonable accuracy, since there is no need to create a fine mesh to resolve separate quantum wells.

### Refractive indices in FDE

Values of refractive indices for InGaAsP layers are taken from [2] for InP lattice matched material. This is a good approximation given the lack of data in literature for specific molar fractions used in this example.

### Frequency range in mqwgain

Since the recombination rates are calculated based on the integral of the spontaneous emission spectra, the frequency range over which the mqwgain solver is run, given by * mqw.lambdamin * and * mqw.lambdamax * , must be wide enough to capture most of the region where spontaneous emission is significant. This ensures that we capture all the spontaneous emission in our estimate of the recombination rates. To reduce the simulation time a modest number of frequency points is set with * mqw.numfreq * and the results are then upsampled using spline interpolation which works very well due to the smoothness of the emission coefficient curves ( * common.frequpsample * and * common.frequpsampletype=1 * ).

### Material properties in mqwgain

Since accurate ternary and quaternary semiconductor material properties are not always available and may be process dependent it is a good practice to tweak and fit some of the properties to match the measurement. In this example specifically, we override the default value of the band gap and include its temperature and carrier density dependence found in [1]. This is implemented in runMQW.lsf file. The material assumed in the master input.lsf file is InGaAsP, which is hard coded in runMQW.lsf in the call to buildmqwmaterial script command. Users can set layer thicknesses, molar fractions, strain, and some other parameters from input.lsf by setting the corresponding options in the mqw struct variable.

### Transverse wave vector range in mqwgain

This wave vector is in the plane of the quantum wells and it is important to use a sufficiently large part of the Brillouin zone to get accurate carrier density in the quantum wells as a sum over the transverse wave vectors. The larger the transverse wave vector the larger the carrier energy in the quantum well and the smaller the probability of occupation of that energy level. With larger temperature higher transverse wave vectors should be included. Typically 10-20% of the Brillouin zone is enough. The value of 20% of the Brillouin zone is hard coded in runMQW.lsf (option simParameters.kt).

### Uncoupled well approximation in mqwgain

If the barriers in the MQW region are thick, for example more than 5 nm, it is a good approximation to consider quantum wells uncoupled by setting * mqw.uncoupled_wells * = true in * input.lsf * . In this case it is sufficient to solve a structure consisting of a single quantum well only, which speeds up the simulation significantly, while producing similar gain and spontaneous emission curves (after appropriate scaling back to the full width case).

### Nonradiative recombination rates in TWLM

SRH, Auger, and spontaneous recombination rates are very important in the overall balance for the laser dynamics. SRH and Auger recombination rates are calculated analytically in script, using well known model equations, based on the average charge density in the MQW region. Since most of the density is in the quantum wells and only a small portion is in the barriers we scale the average input density to get the density in the quantum well and use that in the calculation of the recombination rates. There is a fitting factor * twlm.qw_density_correction * in * input.lsf * where users can adjust the ratio of the density in the wells to the density in the barriers to obtain a better agreement of the threshold current. Spontaneous recombination rate is obtained from mqwgain solver, while the parameters for the SRH and Auger recombination, such as lifetimes and Auger coefficients can be set through the twlm struct variable in * input.lsf * .

### Leakage current in TWLM

The leakage current is that component of the input current that is not injected into the MQW region to contribute to lasing, but leaks through the device. To include this effect the input current is multiplied by a factor smaller than 1 by placing a gain element between the DC source and the TWLM element in INTERCONNECT simulation. The value of this factor is set by * twlm.current_injection_efficiency * in * input.lsf * . User should estimate the value of this parameter. More details on this can be found in taking the model further section.

### Gain and spontaneous emission curve fitting in TWLM

TWLM uses a sophisticated fitting algorithm to include the frequency and density dependence of gain and spontaneous emission coefficients into the time-domain simulation. Users can check the result of fitting and adjust the fitting options if necessary to make sure there are no unphysical peaks in the fitted gain curves that can produce fake lasing thresholds at wrong frequencies. Specifically, the fitting algorithm assumes periodic boundaries, such that the first and the last frequency point have the same gain and spontaneous emission. The rolloff parameter * twlm.fit_rolloff * in * input.lsf * controls the percentage of the total frequency range on both ends and is used to accommodate these periodic boundary conditions. Increasing this value may help obtain a better match at the ends of the frequency range.

## Updating the model with your parameters

Instructions for updating the model based on your device parameters

### Changing the gain material

Modifications must be made to the input.lsf file. An estimate of the Auger recombination coefficients at 0 K as well as their activation energies must be provided. The SRH lifetimes of the electrons and holes in the new well material must be provided. The center frequency and simulation bandwidth of the INTERCONNECT simulation and the center frequency of the FDE simulation must be modified to include the estimated lasing frequency. The frequency range for the mqw solver must be adjusted.

The laser structure must be updated in the eigenmode solver project file (MODE or FEEM) so that the effective and group indices and the confinement factor of the mode to the gain region may be calculated. The group index calculation requires knowledge of the derivative of the effective index. This derivative is calculated by a central difference approximation from two perturbations around the frequency of interest. Since the material dispersion will likely be significant in the group index calculation you must provide script functions that yield the material indices as a functions of frequency and substitute them in * runFDE.lsf * for the function n_InGaAsP__InP_HHI_Eg and then make sure that each material is updated inside the main for loop in runFDE.lsf.

Materials must be created and loaded into the mqwgain solver in the script file * runMQW.lsf * by modifiying a call to to * buildmqwmaterial * script command.

All four script files need to be rerun.

### Changing the thickness of MQW layers

Modifications must be made to the input.lsf file. The quantum well structure must be modified to reflect the new layer thickness. The center frequency and simulation bandwidth of the INTERCONNECT simulation and the center frequency of the FDE simulation may need to be modified to include the estimated lasing frequency. This is due to the quantization effects of energy levels in the quantum wells which depend on the thickness. Similarly, the frequency range for the mqw solver may need to be adjusted.

The laser structure must be updated in the eigenmode solver project file (MODE or FEEM) so that the effective and group indices and the confinement factor of the mode to the gain region may be calculated.

All four script files need to be rerun.

### Changing the number of MQW layers

Modifications must be made to the * input.lsf * file. The quantum well structure must be modified to reflect the new number of layers.

The laser structure must be updated in the eigenmode solver project file (MODE or FEEM) so that the effective and group indices and the confinement factor of the mode to the gain region may be calculated.

All four script files need to be rerun.

### Changing the gain region width

Modifications must be made to the * input.lsf * file. The active region size must be modified to reflect the new width.

The laser structure must be updated in the eigenmode solver project file (MODE or FEEM) so that the effective and group indices and the confinement factor of the mode to the gain region may be calculated.

All four script files need to be rerun.

### Changing the gain region length or the facet reflectivity of the laser cavity

Modifications must be made to the * input.lsf * file. The active region size must be modified to reflect the new length or the left and right facet reflectivities must be modified to reflect the new values.

* runTWLM.lsf * needs to be rerun.

## Taking the model further

Information and tips for users that want to further customize the model

- Carrier-dependent loss due to free-carrier absorption is set to 2200 1/m in this example. Carrier-dependent loss due to free-carrier absorption can be calculated in MODE or FDTD. TWLM can then be updated with calculated result.
- Spontaneous emission coupling factor is set to 0.01 in this example. Spontaneous emission coupling factor can be calculated either using FDTD or based on the numerical aperture of the waveguide. TWLM can then be updated with calculated result.
- To take leakage current into account, laser leakage current at different drive currents can be calculated using 2D/3D CHARGE simulation.
- This Fabry-Perot laser example can be extended to a DFB laser , DBR laser or a Ring Vernier laser following instructions in the linked example.
- Additional laser metrics, not shown here, can be calculated either as a post-processing step on the laser spectrum, such as side mode suppression, or directly visualized from the TWLM solver, such as longitudinal position and time dependent charge, photon density, gain, mode amplitude, and recombination rate profiles.

## Additional resources

Additional documentation, examples and training material

### Related publications

- J. Piprek, P. Abraham, and J.E. Bowers,“Self-Consistent Analysis of High-Temperature Effects on Strained-Layer Multiquantum-Well InGaAsP–InP Lasers”, IEEE Journal Of Quantum Electronics, Vol. 36, No. 3, March 2000, pp. 366-374.
- S. Seifert and P. Runge, “Revised refractive index and absorption of In(1-x)Ga(x)As(y)P(1-y) lattice-matched to InP in transparent and absorption IR-region”, Optical Materials Express ,Vol. 6, No. 2, February 2016, pp. 629-639.

### See also

## Appendix

Additional background information and theory

### Optical Mode Simulation with FEEM

As an alternative option, the laser structure can be created in the Finite Element (FE) IDE for simulation using FEEM. FEEM solver can be used to calculate the optical mode profile and extract the effective index and group index of the fundamental (TE) mode as well as its confinement factor with respect to the gain medium. These are all calculated in the vicinity of the nominal target lasing frequency and are expected to be locally weakly frequency dependent. Once the mode calculation is done, the rest of the steps are identical. Hence, alternatively for step 1:

** Step 1: Optical mode calculation **

- Launch FEEM.
- Open and run the script file input.lsf.This will initialize input data and simulation configuration parameters that are required for all the three steps. This data will be saved in
*Piprek2000OQE.json*which will be called later on by scripts in different steps when needed. - Open and run the script file runFEEM.lsf. This will calculate the effective index as a function of frequency and then the group index using a finite difference approximation.

The effective and group indices, as well as the confinement factor used in the subsequent steps are stored in * wgSweep.neff, wgSweep.ng * , and * wgSweep.confinement_factor * , respectively. These results will also be stored in the file * Piprek2000OQE_eigenmode.json * .

* runMQW.lsf * script in Step 2 can either run in FEEM or INTERCONNECT.