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.

Note: 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 the "Taking the model further" section.

Run and Results

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]

Launch MODE. Open and run the script file input.lsf. This will generate the file Piprek2000OQE.json.
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, wgSweep.ng, and wgSweep.confinement_factor, respectively. These results will also be stored in the file Piprek2000OQE_eigenmode.json.

[Option B, FEEM]

Launch FEEM. Open and run the script file input.lsf. This will generate the file Piprek2000OQE.json.
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, wgSweep.ng, and wgSweep.confinement_factor, respectively. These results will also be stored in the file Piprek2000OQE_eigenmode.json.

[Common Steps]

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.
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).
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.

Results

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

References

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

Taking the Model Further

The workflow can be expanded by using 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.py , which should be run after the TWLM simulation (the last step in the main workflow in the previous post). This python file can be run from INTERCONNECT (or any other product) and it will launch CHARGE solver, create device model and perform the simulation. Python distribution comes pre-installed with our products and the user only needs to ensure the following option is enabled in Help/Python integration status menu:

CHARGE simulation will leverage the results of the FDE, MQW, and TWLM simulations by reading the corresponding result files. runCHARGE.py relies on the following scripts: builddevice.py, mqw_material_build_functions.lsf, mqw_material_database.lsf, 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 files runCHARGE.lsf and builddevice.lsf can be used as an alternative to runCHARGE.py and builddevice.py if the user wants to use MATLAB functions instead of python functions for the least square fit of spontaneous emission coefficients in CHARGE solver from the MQW solver results. This fits the MQW’s spontaneous recombination rate as a function of charge density to a polynomial with the second and third order terms, which is a more accurate model at higher densities where there is a roll off of the recombination rate. There is no difference in results between the .lsf and .py versions of these two files. However, the .lsf files require MATLAB to be installed on the same machine by the user.

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 1V bias:

It is also possible to visualize quantities such as 2D electric field profile:

From the x component of electron current density we see that there is diffusion of carriers out of the active region under the ridge to the region on the right hand side of the ridge, which represents lateral leakage. All electrical results above are visualized at the 1V bias. Additional quantities, not shown above, are available as well (e.g. hole current density).

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.

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

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