Micro-LED

In this example, we simulate a micro-LED. The self-consistently coupled CHARGE and MQW solvers are used to calculate the \(I-V\) response, spontaneous emission power spectrum, and quantum efficiency. The STACK and FDTD are used to characterize and extract the emitted power and radiation pattern.

download example

Overview

Understand the simulation workflow and key results

Similar to OLEDs, micro-LEDs are challenging devices to simulate. They require a large number of simulations and some post-processing is needed to analyze the devices' performance.

Here we will use the CHARGE/MQW self-consistently coupled solvers to simulate the IV curve, spontaneous emission power spectrum, and quantum efficiency. We will also use the STACK optical solver to characterize the emission properties of a planar device, and FDTD to simulate a 3D micro-LED.

Step 1: Simulate the micro-LED using Coupled CHARGE/MQW Solvers

In this step we simulate the electronic and optoelectronic figures of merit, such as the IV curve, spontaneous emission power spectrum, and internal quantum efficiency.

CHARGE/MQW coupled solver ensures self-consistent coupling between bulk semiconductor domains with solutions of the drift-diffusion and Poisson’s equations solved by CHARGE and quantum well active domain with solutions of the Schrödinger equation solved by MQW using k.p method..

Step 2: Characterize the 1D stack using STACK Optical Solver

The incoherent, unpolarized emission from the active layer of the planar LED is obtained using the STACK optical solver.

This approach can be used to quickly characterize and optimize the layers' thickness and materials before simulating the full 3D device with FDTD.

Step 3: Simulate the cylindrical micro-LED using 3D FDTD

We use a series of FDTD simulations to extract the incoherent, unpolarized emission of the micro-LED and determine the light extraction efficiency.

Run and Result

Instructions for running the model and discussion of key results

Step 1: Simulate the micro-LED using CHARGE/MQW Coupled Solver

1. Open file [[Horng2018.ldev]] and run CHARGE solver to simulate the micro-LED.
2. Open and run script [[Horng2018_plot.lsf]] to postprocess and plot results.

In the first step, the project file is already set up with all the material, geometry, and solver parameters that correspond to the device from reference [1] (Horng, 2018). Running the CHARGE solver with the MQW domain option set in CHARGE edit window will ensure running CHARGE and MQW in the self-consistently coupled mode. MQW options can be set from the CHARGE solver edit window, and they are similar to the standalone MQW options, with unused options in the coupled mode being greyed out.

Please note that this is a time consuming simulation, due to the existence of 10 quantum wells that are solved quantum-mechanically multiple times to achieve convergence of total density, potential, and recombination rates. To speed up the simulation as much as possible we use a narrow simulation domain, effectively making the simulation 1D. This is a good approximation, but the current spreading effects are not included due to the narrow simulation domain.

In the second step, we post-process the results and plot the main figures of merit. In order from left to right and top to bottom they are in the figures below: IV curve, spontaneous emission total power, spontaneous emission power spectrum, and the quantum efficiency. We also plot measured results from reference [1] for comparison. The simulated IV curve has a similar turn on point as the measured curve, but we assumed a 2 ohm series contact resistance to get a better match for the slope. The simulated spontaneous emission power spectrum was plotted at the highest simulated voltage of 2.5 V.

Step 2: Characterize the 1D stack using STACK Optical Solver

1. Open the simulation file [[micro_LED_cylindrical.fsp]] in FDTD
2. Open and run the script file [[STACK_micro_LED.lsf]]

In this script, we use the STACK optical solver to analyze the planar emission of the device. The script constructs the layers and extracts material properties from the material database.

The script uses the stackpurcell command to calculate the Purcell factor as well as the power density per steradian emitted into the sapphire as a function of angle and wavelength. It also integrates the power density to calculate the total upward power (normalized to the power the dipole would radiate in a homogeneous medium). These results will be easy to compare with the 3D FDTD simulation.

Additionally, we use stackdipole to get the luminance and radiance as a function of angle for a given dipole spectrum, as well as the X, Y, Z tristimulus values.

Step 3: Simulate the cylindrical micro-LED using 3D FDTD

1. [[step]]Open and run the simulation file micro_LED_cylindrical.fsp
2. Open and run the script file plot_single_dipole_results.lsf to plot the following figures

We start by running a first simulation at low accuracy (mesh accuracy set to 1), with a source position set to 0, oriented along the x-axis. We can then visualize the field intensity in the cross-section from the monitor y-normal (shown here in log scale), or from the monitor farfied::z2 (shown in linear scale).

The “farfield” analysis script will calculate the total power flowing out of the box, as well as the normalized Poynting vector in the far-field (at a distance of 1m). This is normalized so the integral of the Poynting vector is equal to the total power propagating in the entire upper half-space, normalized to the power that would be emitted by the dipole in a homogeneous medium of the same refractive index as the MQW layer. In addition, the analysis group calculates P_vs_theta, the average value of the Poynting vector over all azimuthal angles, and, therefore, represents the Poynting vector as a function of polar angle for a ring of incoherent dipole emitters at the same radial position from the center of the device.

Lastly, it calculates the Purcell factor on the same wavelength grid as the other results (the Purcell factor can be easily calculated on any wavelength/frequency grid from the dipole source).

These results can be then used to perform further calculations. For example, we calculate the total extraction efficiency into a cone of 34.4 degree (extraction from sapphire to air):

3. [[step]]Run the parameter sweep "position"
4. Open and run the script file dipole_position_analysis.lsf

The nested sweep stores all the results of the “farfield” analysis group for various positions and polarizations of the dipole. To save the memory, we can remove the far-field result which consumes most memory.

The sweep is set to calculate 11 positions for each of the 3 dipole polarizations, for a total of 33 simulations. To reduce the simulation time, run this sweep with a mesh accuracy set to 1, and set the property “project all wavelengths” to 1 (true).

The script dipole_position_analysis.lsf will take the value of P_vs_theta and take the average of the 3 dipole orientations to calculate the result for unpolarized dipole emission (see the " Important model settings " ). Then, it will calculate weights for the contribution of each dipole position, assuming a uniform dipole density per unit area. Each dipole is therefore multiplied by a weight that is proportional to the area of the ring it represents, as shown in the following figure:

These weights can be easily adjusted if more information is known about the current density in each region.

The script calculates the power density vs theta at a single wavelength near 625nm or spectrally integrated over a dipole spectrum about 10nm wide near the same wavelength.

It also calculates the total power emitted in a 34.4 degree cone where it can be emitted to air.

Important model settings

Description of important objects and settings used in this model

CHARGE electrical material properties

The micro-LED device in this example is made of III-V ternary and quaternary compounds. The electrical material database stores default material properties and their interpolation for arbitrary compound composition. However, it may be necessary to check important options and tweak their values in case the default interpolation is not good enough. Especially important properties to confirm are band gaps, mobilities, and SRH (trap) and Auger recombination rates. Specifically, recombination rates are dependent on the process and cannot be accurately obtained from interpolation.

In this example, the main radiative recombination rate is located in the MQW active layer and is calculated quantum-mechanically. However, in the material properties there is also the bulk radiative recombination parameter. This parameter controls radiative recombination in bulk regions. In this example, we set the bulk radiative recombination parameter to 0 in the SCL and AlInP layers to prevent spurious emission effects.

CHARGE doping

The doping profile is very important, especially for the IV curve and device turn on. In this example, the doping profile was taken from reference [1].

MQW band offsets

The band offsets are very important both for energy sub-bands and transition in quantum wells as well as for current thermionic leakage. Default band offsets in the electrical material database may not always be accurate for every material and composition. It is advisable to check band offsets and tweak them using custom models available in the literature. Band offsets can be set by setting the work function in the electrical material database as explained in this article Understanding band offsets in heterostructures – Ansys Optics . They are applied to both CHARGE and MQW simulation.

MQW parameters

The rest of MQW parameters should be set in the CHARGE solver edit window in a similar way to how they are set in standalone MQW gain simulations in our laser simulation examples. Specifically, the spatial grid number of points, transverse wave vector number of points, and number of frequency points affect the MQW simulation time. These options should be set to the lowest acceptable values in order to speed up CHARGE/MQW coupled mode simulations.

FDTD "farfield" analysis group

The " farfield " analysis group has a setup script to create a box of monitors of the desired size, which can be positioned over the LED structure.

Since it is on a GaAs substrate, we can choose not to include the monitor in the z min boundary (property "include_z_min_boundary" set to 0).

The following analysis properties can be modified by the user:

• theta min, theta max, N theta : min, max, and number of points for polar angle theta.
• N phi per half : number of values of phi between 0 and 180 degrees. The Poynting vector from 180 to 360 degrees is calculated assuming symmetry.
• project all wavelengths : can be 0 (false) or 1 (true). If true, all wavelengths will be projected which can take 10s of seconds or more. If false, the single wavelength closest to the " target wavelength " will be chosen for projection.
• near field points per wavelength : the near fields can be down-sampled for far-field projection, which substantially speeds up the projections. Values of 3 or 4 are typically sufficient. Increasing this will result in longer far-field projection times without improving accuracy.

FDTD mesh accuracy

In the above steps, we use a mesh accuracy of 1 to speed up the simulations. For better accuracy, we can set up the mesh accuracy to 2 or 3 (See the " Convergence " ).

FDTD dipole orientation

The emission of such micro-LED is incoherent and unpolarized. While FDTD is fundamentally a coherent simulation method where sources have a well-defined polarization, we can extract incoherent, unpolarized results by adding incoherently the fields from simulations with orthogonal polarizations. This means we have to run 3 simulations for each dipole position, 1 for each possible dipole orientation, x, y, and z (see Incoherent, unpolarized dipole ).

The dipole orientation is set in the "model" object. This group is also parametrized for the dipole position along the x-axis. It will automatically adjust the symmetric and anti-symmetric boundaries of the simulation, assuming the dipole is located on the x-axis, at \(y=0\) to reduce the simulation volume where possible. The properties of " model " are shown below.

The model setup script will force the x, y position and orientation of the dipole source, and set the x min, y min boundary conditions of the FDTD region. These changes will overwrite any direct changes made to those properties.

Updating the model

Instructions for updating the model based on your device parameters

FDTD

The micro-LED structure is parametrized in a structure group that allows us to easily change the layer thicknesses and materials.

The dipole is centered in the middle of the MQW layer. This position can be changed by modifying the value of "z" in the source properties.

Taking the model further

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

FDTD materials

The optical materials used in this example are taken from the reference paper. Other materials can be used provided you know their optical properties over the bandwidth of interest. See the appendix for more information about the materials used and how they were added to the FDTD simulation.

FDTD geometry

The device studied here is a cylindrical micro-LED. Other geometries can be used, but note we use the symmetries extensively to reduce the simulation size as well as the number of dipole positions. This should be taken into consideration when updating the geometry.

FDTD dipole spectrum

FDTD will normalize the fields recorded by frequency monitors to the source spectrum, so these monitors will return the impulse response of this device. It is then easy to calculate the response to an arbitrary power spectrum. We simply multiply the impulse response of the system by the new power spectrum:

$$ f_{spectrum} = Power_{density}(\lambda) \cdot f_{impulse}(\lambda) $$

Where \(f_{impulse}\) is the impulse response and \(Power_{density}\) the new power spectrum.

This can be performed as a post-processing step. Therefore, you can calculate the response to any spectrum without having to re-run the FDTD simulations.

Export to ray tracing

Designing OLED/LED devices requires a combination of nanoscale and macroscale optics. The individual pixels often have sub-wavelength features such as thin dispersive layers, scattering structures, that require an electromagnetic field solver, while the emission from the macroscopic device requires a ray-based tool.

For these applications, one can calculate the angular distribution of a pixel with the STACK optical solver or FDTD and then load the result into a ray-tracing tool in the form of ray sets.

FDTD supports the export of the angular distribution to ray sets for Breault Research ASAP, and ZEMAX OpticStudio (requires an additional license).

Convergence Test

Tips for ensuring that your model is giving accurate results

The script [[ dipole_position_analysis.lsf]] makes it possible to use less than 11 dipole positions originally simulated, to understand if fewer positions could be used. In the following figures, we compare 11 dipole positions, with 6 and 3 positions only.

This comparison shows that 6 dipole positions may be sufficient for many analyses, while 3 is not.

We can also run the parameter sweep with a FDTD mesh accuracy of 2 (targeting \(dx=\lambda/10\) instead of \(dx=\lambda/6\)) in all media. Mesh accuracy of 1 is sufficient for rapid simulation and optimization, but errors of 10 to 20% are typical due to the extremely coarse FDTD mesh size. The change from mesh accuracy 2 to 3 (\(dx=\lambda/14\)) is typically only 1-2%.

Additional resource

Additional documentation, examples, and training material

Related publications

1. R. Horng, H. Chien, K. Chen, W. Tseng, Y. Tsai and F. Tarntair, “Development and Fabrication of AlGaInP-Based Flip-Chip Micro-LEDs,” in IEEE Journal of the Electron Devices Society, vol. 6, pp. 475-479, 2018 (doi: 10.1109/JEDS.2018.2823981).

See also

• Far-field projections in FDTD overview
• OLED (2D)

Related Ansys Innovation Courses

Appendix

Additional background information and theory

2D Simulations

A 2D simulation of micro-LED can be useful if you want to quickly check its behavior at the initial design stage.

Simulation setup

The micro-LED is parameterized in a structure group which allows the material thicknesses and materials to be easily changed. For 2D simulations, the plane of the simulation is the x-y plane, therefore the orientation of the structure is such that the normal to the surface is the y-axis. The dipole is centered in the middle of the MQW layer. The entire model group is parameterized for the dipole position (along the x-axis) and axis (1, 2 or 3 for x, y and z polarized dipoles). The structure is assumed to be symmetric with respect to the x=0 plane. Therefore the "model" group will also automatically adjust the symmetric and anti-symmetric boundaries of the x-min boundary when the dipole is placed at x=0 to reduce the simulation volume. The properties of the "model! are shown below. The model setup script will force the x, y position and orientation of the dipole source, and set the x min, y min boundary conditions of the FDTD region. These changes will overwrite any direct changes made to those properties.

Single simulation analysis
Run the simulation with the source position set to 0, and the axis set to 1 (x-axis). This will take a few seconds on most machines. After running, we can plot some of the field cross-section, like z-normal, shown on a log scale:

The farfield analysis group is designed specifically for this 2D simulation. It has a setup script to create a box of monitors of desired size, which can be positioned over the OLED structure. Since the structure is on a GaAs substrate, we can choose not to include_y_min_boundary. The farfield analysis group has an analysis script that will calculate the total power flowing out of the box, as well as the normalized Poynting vector in the far field at a distance of 1m. This is normalized to the power that would be emitted by the dipole in a homogeneous medium of the same refractive index as the MQW layer. In addition, the farfield group assumes that the structure is symmetric and therefore calculates the Poynting vector as if there were a second dipole at the -x position of the actual source. This reduces the number of simulations by a factor of 2 when sweeping the dipole position. The farfield monitor will also integrate the power in the far field for comparison to the near field. Since the near field is not a closed box (because we are not including the y-min boundary), and because there is power radiating to the lower half space, these are not expected to agree. The farfield monitor also calculates the integral of the power in the far field weighted by \(sin(\theta)\) to approximate what might be observed in 3D if the 2D results are assumed to be a \(P\) vs \(\theta\) for a cylindrically symmetric 3D structure.
Finally, the farfield group calculates the Purcell factor on the same wavelength grid as the other results. The Purcell factor can be calculated easily on any wavelength/frequency grid from the dipole source.

The following analysis properties can be modified by the user:

• theta min, theta max, N theta : min, max and number of points for polar angle theta. This is normally 0 to 90 degrees because of the assumption of symmetry noted above.
• near field points per wavelength: The near fields can be downsampled for far field projection, which substantially speeds up the projection. Values of 3 or 4 are typically sufficient. Increasing this will result in substantially longer far field projection times without improving accuracy, although the projection times are typically negligible in 2D.

P vs. theta/wavelength:

Transmission vs. wavelength:

The Purcell factor (which can also be obtained to higher resolution directly from the dipole source or the script command dipolepower ):

These results can be used to perform various calculations. For example, if we want to know the total extraction efficiency into a cone of 34.4 degrees (extraction from sapphire to air), we can use the script file [[power_in_angular_cone_2D.lsf]] to generate the following result, where we weighted by \(sin(\theta)\) to better approximate a 3D result.

Parameter sweep of dipole position and polarization
The Optimization and Sweeps window has a nested parameter sweep that can sweep over dipole positions and polarizations. For each dipole position and polarization, it stores all the results of the farfield analysis group. We calculate all dipole polarizations but show later how only in-plane dipoles can be included in the analysis.
Initially, the sweep is set up to calculate 11 dipole positions for each of the 3 dipole polarizations, for a total of 33 simulations. In 2D this runs in a few minutes on most machines.
After completing the 33 simulations, the concatenated results will be stored in the “position” sweep object. If you save the FDTD file, this sweep information will be saved to the file and can be recovered even if you close the file.
Open the script [[dipole_position_analysis_2D.lsf]]. It will take the value of P vs theta, and take the average of the 2 or 3 dipole orientations to calculate the result for unpolarized dipole emission (see the Incoherent Unpolarized Dipole page for more details on this calculation). To include all 3 dipole orientations, comment line 8. It will then calculate weights for the contribution of each dipole position assuming a uniform dipole density per unit area and a cylindrically symmetric structure. This is an approximation because we have only 2D simulation results, but it can more closely match the 3D results for symmetric structures. Each dipole is therefore multiplied by a weight that is proportional to the area of the ring it represents, as shown in the following figure.

These weights can easily be adjusted if more information is known about the current density in each region. In addition, when integrating the Poynting vector in the far field, we weight the angles by sin() to better approximate a 3D result (as if our 2D results for \(P\) vs \(\theta\) are simply a cross-section of a cylindrically symmetric structure). The power density vs theta at a single wavelength near 625nm or spectrally weighted over about 10nm near 625nm is calculated:

The total power emitted in a 34.4-degree cone where it can be emitted to air is shown below.

Convergence
The script makes it possible to use less than the 11 dipole positions originally simulated, to understand if fewer dipole positions could be used. We will study this more carefully in the 3D example when we are trying to reduce the number of simulations as much as possible.

Materials

For materials, we use:

• (AlxGa1-x)0.5In0.5P latticed matched to GaAs from Kato, Hirokazu & Adachi, Sadao & Nakanishi, Hiroshi & Ohtsuka, Kouji, Optical Properties of (AlxGa1-x)0.5In0.5P Quaternary Alloys. Japanese Journal of Applied Physics. 33. 10.1143/JJAP.33.186 (1994).
  We modified the imaginary part of the permittivity with a Heaviside-like function to better model the transition to transparency near 600nm. The script file [[ micro_LED_AlGaInP_mdf.lsf]] will create the necessary models and add them to the current material database. It will remove any pre-existing materials of the same name. It also reproduces Figure 10 of Kato et al:

• GaAp from https://refractiveindex.info/?shelf=main&book=GaP&page=Aspnes (yaml file included).
• Sapphire (Sellmeier model) from https://refractiveindex.info/?shelf=3d&book=crystals&page=sapphire . Note that in FDTD simulations this will use the refractive index value at the center frequency of the simulation and will be non-dispersive.
• ITO from https://refractiveindex.info/?shelf=other&book=In2O3-SnO2&page=Konig .

The materials used have been saved in micro_LED_materials.mdf . The default FDTD fit settings on some of the materials have been modified from their defaults to improve fitting over 550-700nm for FDTD simulation:

• ITO and GaP use a fixed fit range of 400-700nm
• MQW material uses a maximum number of coefficients of 1. This forces it to use a simple lossless dielectric model. The lossless nature is useful for calculating the Purcell factor in FDTD.

Far-field projections

The farfield projection is typically performed from a single surface which is fast because it can use ffts. However, the fields must decay to 0 at the edges of the surface. While filtering of the near field can be used to reduce these effects it still can require larger simulation volumes. If we look at the near field at "farfield::z2", we see that the field is not zero near the edges of the monitor.

A comparison of projection methods below shows that the filter of the near field (farfieldfiler > 0) cannot be used to recover accurate results. The method of projecting from the 5 sides of the box gives the best results, as shown in the bottom right.

• On the left, fft based projection from "farfield::z2" only, with "farfieldfilter" = 0; on the right, fft based projection from "farfield::z2" only, with "farfieldfilter"=1

• On the left, far-field exact projection from only "farfield::z2" (gives same results as fft based projection); on the right, current projection from 5 sides of the box.

In principle, the far-field projection needs to be from a closed surface in a homogeneous medium. In this case, due to the substrate, we have excluded the lower z surface. As a result, the power integrated over the open box in the near field and the power integrated over the upper half-space in the far-field are not expected to perfectly agree. Inaccuracies in the radiated power will predominately affect very steep emission angles, far from normal. Here is a comparison of the two integration methods: