Simulation of LIV Curves for InGaAsP/InP MQW Ridge Laser


A key performance characteristic of semiconductor lasers is the L-I curve, which is a plot of the output light from the laser vs the current injected. In this example we demonstrate a workflow for simulating the L-I curve of an InGaAsP-InP multiple-quantum well (MQW) ridge laser presented in [1].

The input data and simulation configuration parameters are entered into a master input file (Lumerical script file). The laser structure is created in the FD IDE for simulation with MODE and as an alternative option in the FE IDE for simulation using FEEM. A test bench is set up in INTERCONNECT using the TWLM element for the measurement of the LI curves. The workflow consists of running three additional script files which perform the following three steps. Please note that both the MQW gain solver and TWLM used in this example require separate licenses, and that the minimum version of the Lumerical suite is 2019b-R2.

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 or the FEEM solver.

Step 2: Gain Medium Simulation

Using the MQW Gain Solver running in the FEIDE, a 4x4 k dot 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.

Step 3: 1D Traveling-Wave Laser Simulation

Using the TWLM element running in INTERCONNECT a 1D laser simulation is performed using a sweep over different drive currents. The optical power emitted may then plotted as a function of drive current to generate the L-I curve. This is done by right clicking on the sweep, selecting power, and then visualize.


For an example of how to expand this LI simulation workflow with our CHARGE drift-diffusion solver to simulate voltage and leakage current and obtain a full LIV curve refer to one of the following posts in this topic.

Run and Results


input.lsf (6.3 KB)
mode.lms (289.6 KB) runFDE.lsf (8.0 KB)
feem.ldev (5.3 MB) runFEEM.lsf (10.0 KB)
runMQW.lsf (11.2 KB)
twlm.icp (263.9 KB)runTWLM.lsf (14.6 KB)

Steps to Follow

The initial steps involve calculating and analyzing the optical mode of the MQW structure. This can be done by using either the FDE solver in MODE [Option A] or the FEEM solver [Option B]. Once the mode calculation is done, the rest of the steps are identical.

[Option A, FDE]

  1. Launch MODE. Open and run the script file input.lsf. This will generate the file Piprek2000OQE.json.

  2. While still in MODE, 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 different approximation. The effective and group indices, as well as the confinement factor used in the subsequent steps are stored in wgSweep.neff,, and wgSweep.confinement_factor, respectively. These results will also be stored in the file Piprek2000OQE_eigenmode.json.

[Option B, FEEM]

  1. Launch FEEM. Open and run the script file input.lsf. This will generate the file Piprek2000OQE.json.

  2. While still in FEEM, 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 different approximation. The effective and group indices, as well as the confinement factor used in the subsequent steps are stored in wgSweep.neff,, and wgSweep.confinement_factor, respectively. These results will also be stored in the file Piprek2000OQE_eigenmode.json.

[Common Steps]

  1. While still in MODE or FEEM, open and launch the script file runMQW.lsf. This will first load 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 for the MQW gain region and run the script-based MQW solver mqwgain 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.

  2. You can now optionally close MODE or FEEM. Launch INTERCONNECT. Open and run the script file runTWLM.lsf. This will load 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 the gain, spontaneous emission spectrum, and 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. When the sweep completes, the optical power emitted may then be plotted as a function of drive current to generate the L-I curve. This is done by right clicking on the sweep, selecting power, and then visualize. This is the L-I curve for the laser operating at 293K (20C).

  3. To produce the L-I curve at 333K (60C), open the input.lsf and set the parameter twlm.temperature to 333. Then find the parameter twlm.li_currents. There will be a line immediately below it with a list of currents more appropriate for the producing the L-I curve at 333K. Use those currents. Run the input.lsf file, then run the runTWLM.lsf file. When the INTERCONNECT simulation has completed you can visualize the L-I from the sweep as in Step 4.


Figure 1 depicts the simulated L-I curves at two different temperatures. The threshold currents and slopes from the present Lumerical simulation as well as those from Reference 1 are presented in Table 1.

Figure 1


Table 1

Lumerical Reference 1
Threshold Current {A] 273K 0.120 0.114
Slope [W/A] 273K 0.27 0.24
Threshold Current [A] 333K 0.250 0.247
Slope [W/A] 333k 0.23 0.19


[1] “Self-Consistent Analysis of High-Temperature Effects on Strained-Layer Multiquantum-Well InGaAsP–InP Lasers”,J. Piprek, P. Abraham, and J.E. Bowers, IEEE Journal Of Quantum Electronics, Vol. 36, No. 3, March 2000, pp. 366-374.

[2] 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.

Important Model Settings

Since the recombination rates are calculated based on the integral of the spontaneous emission spectra, the frequency bandwidth over which the mqwgain solver is run must be wide enough such that the spontaneous emission curves become insignificant near the mqwgain simulation band edges. This ensures that we capture all the spontaneous emission in our estimate of the recombination rates.

Updating the Model with Your Parameters

You may adapt this example to use different materials and a different layer structure for the layers. This will require the laser structure to be updated in the eigenmode solver (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 (or runFEEM.lsf) for the function n_InGaAsP__InP_HHI_Eg and then make sure that each material is updated inside the for loop in runFDE.lsf at line 294 (or runFEEM.lsf at line 361).

Modifications must also be made to the input.lsf file. The quantum well structure must be modified to reflect the new quantum well structure, e.g., the number and width of the of the wells must be modified An estimate of the Auger recombination coefficients at 0 K as well their activation energies must be provided. The Shockey-Read-Hall lifetimes of the electrons and holes in the new well material must be provided. The center frequency and simulation bandwidth of the INTERCONNECT simulation must be modified to include the estimated lasing frequency. Materials must be created and loaded into the mqwgain solver in the script file runMQW.lsf

Taking the Model Further

Estimating the spontaneous emission coupling either using FDTD or a based on the numerical aperture of the waveguide.

Calculating the carrier-dependent loss due to free-carrier absorption in FDTD. The carrier dependent loss can then be subtracted from the gain before creating the files containing the gain spectrum.

Calculating the leakage current at different drive currents using CHARGE.

Additional Resources





I ran the script file input.lsf and mode solutions is showing an error
input (1).lsf line 84: jsonsave is not a valid function or variable name
please help


Do you have a recent enough version of our products ? Maybe your version does not have the jsonsave, jsonload script commands. Try running the following the commands in the script prompt :


The previous workflow for calculating LI laser curves, where I is the total recombination current in the active region, can be expanded by using our CHARGE drift-diffusion solver to include electrical contacts and other domains beside the active region, such as ridge, cladding, and separate confinement layers. CHARGE simulation enables calculating the total current, including leakage current, as well as the voltage in order to obtain full LIV curves. For this purpose we will use the new stimulated recombination model for semiconductor materials in CHARGE solver, which can be set by using the stimulated recombination rate results calculated by the TWLM solver.

The main simulation file is runCHARGE.lsf (6.6 KB), which should be run after the TWLM simulation. CHARGE simulation will leverage the results of the FDE, MQW, and TWLM simulations by reading the corresponding result files. runCHARGE.lsf relies on the following scripts: builddevice.lsf (22.5 KB) mqw_material_build_functions.lsf (10.5 KB) mqw_material_database.lsf (1.7 KB), which are used to automatically create device model (e.g. quaternary material properties, device geometry, solver options) in script and should be stored in the same folder as the main files.

NOTE: Script file runCHARGE.lsf calls certain MATLAB functions to fit the MQW’s spontaneous recombination rate as a function of charge density to a cubic polynomial required by the spontaneous recombination model in CHARGE. If you would prefer to avoid calling MATLAB and make Lumerical script self-contained let us know by posting in this topic.

After running the CHARGE simulation with runCHARGE.lsf the following results will be obtained. The LIV curve
Total current and different components of current

I*dV/dI curve showing the threshold point where the curve has a sharp drop
The users can also visualize the stimulated recombination rate such as the plot below for above threshold bias of 1V

where we can see that the stimulated recombination is nonzero only in the thin active region below the ridge as expected.

Since this is a 2D CHARGE simulation, in order to reduce the simulation time we reduced the total width of the ridge and substrate to 8 and 16 microns, respectively, and simulated only half of the structure by exploiting symmetry. The total current for the actual full width can be found by scaling the calculated current to the total width, exploiting the relative uniformity of current density in the ridge away from the edges, and by multiplying by factor 2 to account for simulating only half of the structure.

Since the active region is much thinner then the rest of domains, we represented the 6 actual quantum wells with a single quantum well of the total thickness equal to the thickness of combined quantum wells. In this way the finite element mesh will be easier to generate and will have less vertices, reducing the simulation time, while still being electrically equivalent to the original case because the charge density in the active region will be unchanged due to the unchanged quantum well volume.



I am trying to run the code runCHARGE.lsf after all the steps and generating all the files required. I also integrated MATLAB into my system.

However, I cannot succesfully run the code since its showing an error at line 184 for buildDEVICE telling:
struct does not contain a field named strainedBandgaps
Could you please address this issue?