In this example we demonstrate a workflow for simulating the LI curve of an InGaAsPInP multiplequantum well (MQW) edge emitting laser (EEL). Calculated LI curves and lasing frequencies at different temperatures are compared to results reported in literature.
MODE MQW INTERCONNECT Laser Library
Minimum product version: 2019
Updating the model with your parameters
Understand the simulation workflow and key results
Laser library and MQW module are used to model a FabryPerot InGaAsPInP MQW ridge laser using scripts. LI 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 IIIV semiconductors with zincblende crystal structure. The simulations are isothermal, suitable for lower power or pulsed lasers without pronounced selfheating 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.
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.
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.
Using the TWLM element running in INTERCONNECT a timedomain 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 LI curve. The spectrum can be plotted at different times.
Instructions for running the model and discussion of key results
1.Launch MODE.
2.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.
3.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_eigenmode.json.
1.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).
1.Launch INTERCONNECT.
2.Open and run the script file runTWLM.lsf.
3.Right clicking on the sweep, selecting power, and then visualize. This is the LI curve for the laser operating at 293K (20C).
4.Right clicking on the sweep, selecting spectrum, and then visualize. For a desired current, spectrum of the laser can be plotted at 293K (20C).
5.Open the input.lsf and set the parameter twlm.temperature to 333 to produce the LI curve at 333K (60C).
6.Find the parameter twlm.li_currents and uncomment the list of currents more appropriate for the producing the LI curve at 333K.
7.Run the input.lsf file again.
8.Run the runTWLM.lsf file.
9.Right clicking on the sweep, selecting power, and then visualize. This is the LI curve for the laser operating at 333K (60C).
10.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 nonradiative 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. 34), 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 (multithreading) may not always produce the desired speedup.
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 LI 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 
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.
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.
Values of refractive indices for InGaAsP layers are taken from ref. [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.
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).
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 ref. [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. The 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.
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 1020% of the Brillouin zone is enough. The value of 20% of the Brillouin zone is hard coded in runMQW.lsf (option simParameters.kt).
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).
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 the user 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 an Auger coefficients can be set through the twlm struct variable in input.lsf.
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. The user should estimate the value of this parameter. More details on this can be found in taking the model further section.
TWLM uses a sophisticated fitting algorithm to include the frequency and density dependence of gain and spontaneous emission coefficients into the timedomain simulation. The 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.
Instructions for updating the model based on your device parameters
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.
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.
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.
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.
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.
Information and tips for users that want to further customize the model
▪Carrierdependent loss due to freecarrier absorption is set to 2200 1/m in this example. Carrierdependent loss due to freecarrier 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 FabryPerot 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 postprocessing 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 documentation, examples and training material
[1]J. Piprek, P. Abraham, and J.E. Bowers,“SelfConsistent Analysis of HighTemperature Effects on StrainedLayer MultiquantumWell InGaAsP–InP Lasers”, IEEE Journal Of Quantum Electronics, Vol. 36, No. 3, March 2000, pp. 366374.
[2] S. Seifert and P. Runge, “Revised refractive index and absorption of In(1x)Ga(x)As(y)P(1y) latticematched to InP in transparent and absorption IRregion”, Optical Materials Express ,Vol. 6, No. 2, February 2016, pp. 629639.
Useful links to documentation
Alternative approach to calculate optical mode profile in Step1.
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
1.Launch FEEM.
2.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.
3.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.