This example performs 3D simulation of a GaN Nanowire (NW) LED with an embedded InGaN quantum disk (QD).

We first perform a strain simulation, to get deformation potentials and piezoelectric polarization, than we apply drift-diffusion model with an increasing bias to the contacts, until the nanocolumn diode is brought in conduction regime. Then quantum efa calculations are performed to get the electron and hole states in the QD.

The device structure is defined in the geometry .geo file and is the following:
1 µm GaN NW, with a radius R = 25 nm
5 nm thick InGaN QD with 10% In content
10 nm GaN barriers on top and below the QD

The NW LED is a nip structure oriented along (000-1) cristallographic direction

While InGaN QD active region is kept intrinsic, p-GaN region is doped with Mg acceptor doping, 7x10¹⁸ cm^-3 and n-GaN is doped with Si donor doping, 3x10¹⁹ cm^-3

In order to execute correctly this example you should have the following files in the same working directory:
InGaN_NW.tib: input file for TiberCAD
InGaN_NW.msh: mesh file produced by the GMSH script InGaN_NW.geo

DEVICE STRUCTURE

In the following, some features of the input file will be described ; for further details you can refer to the program reference manual.
In the Device section, the structure of the nanocolumn is described: the InGaN quantum disk region, the GaN barriers and the two high doped GaN contact regions.
First, the crystal directions are defined:

structure = wz
x-growth-direction = (1,1,-2,0)
y-growth-direction = (1,-1,0,0)
z-growth-direction = (0,0,0,-1)

Note that the crystal structure is defined as wurtzite

structure = wz

The crystal orientation is defined only by the 2 vectors

x-growth-direction
y-growth-direction

while

z-growth-direction

is derived internally from x and y vectors. This means that the orientation along z-axis must be checked by user to be consistent with a right-hand convention with respect to x and y. In this case the direction along z is opposite to this, meaning that we have set up a (000-1) orientation.

When in the modeling script file InGaN_NW.geo we have defined

Physical Volume(qdisk) 

we have created in the mesh a physical region with the name

qdisk

we can refer to this GMSH physical region (which is a volume, in this 3D example), by writing

Region qdisk

In the same way we can refer to the GaN volume regions, for example for the one defined with:

Physical Volume(p_Gan)

We write

Region p_GaN 
{
material = GaN 
Doping
{ density= 1e18
type = acceptor
level = 0.17 
}
}

Energy level of the dopant (Mg in this p-doping case) is defined by

level = 0.17

(it determines incomplete dopant ionization at room temperature).

SIMULATION MODULES

1. Strain

We want strain calculation, so we define a simulation

name = strain

belonging to the model elasticity, to be associated to the whole device

regions = all

a

lattice_mismatch model

is defined, with GaN as reference material:

Module elasticity{ 
 name = strain regions = all
} 
...................... Physics 
{ body_force lattice_mismatch 
{
reference_material = GaN 
x-growth-direction = (1,1,-2,0)
y-growth-direction = (1,-1,0,0)
z-growth-direction = (0,0,0,-1)
}
} Contact anode {type=clamp}

2. Drift-diffusion

As for drift-diffusion, as usual we define a simulation

name = dd

belonging to the model driftdiffusion and associated to the whole device

regions = all

A srh recombination model is implemented and also one for radiative (direct) recombination, in this way:

recombination srh 
{}

recombination direct
{}

Field dependent mobility model is chosen for both holes and electrons.

The Boundary Regions for drift-diffusion are the two contact regions, defined by the two boundary surfaces anode and cathode

In Physics section it is in general possible to choose a model with which to obtain band properties (see User Manual for details).
Band parameter models may be defined through the block density_of_states.
This can be done inside:

1. a single block named band_properties containing parameters for both conduction and valence band
2. separate blocks named conduction_band and valence_band to control each band independently

Since in this example we want the strained implementation of drift-diffusion model, we have to define a band_properties block where we choose a density_of_states model for the calculation of bulk band parameters based on kp, for both conduction and valence band.

band_properties 
{
 density_of_states bulk_kp
 {
 strain_simulation = strain
 }
}

This model, defined by

density_of_states bulk_kp

calculates the band edge energies and masses from bulk kp theory, including Pikus-Bir strain corrections. To include strain corrections, the keyword strain_simulation has to be used, providing the name of the module instance which calculates strain. In the current version it is not possible to provide kp parameters from the input file.

Important
In the previous TiberCAD releases the bulk_kp model for the DOS was the default model.
Now the default model is the Bulk  model, which describes the 3D DOS of a parabolic band semiconductor. Band edge energies and the DOS mass are taken from the database. This is particularly suited for Si.

However, when material different from Silicon are present in the device, or when strain has to be included in calculations, as it is the case of this example, it is necessary to define the bulk_kp model for density_of_states.

Moreover, we need to include  the piezoelectric polarization (besides pyro or spontaneous  polarization) arising from strain in the wurtzite nitrides materials of this nanocolumn.  In fact it will enter in Poisson equation and heavily modify the band profiles. It is calculated with:

polarization piezo { strain_simulation = strain }
polarization pyro {  }

Note that if

strain_simulation = strain

is not defined, strain results, even if calculated, will NOT be used in drift-diffusion model.

3. Quantum EFA

For the quantum calculations, we define two efaschroedinger simulations:

quantum_el

for electrons and

quantum_hl

for holes

For electrons we write:

Module efaschroedinger
{
 name = quantum_el
 regions = quantum 
 plot = (ProbabilityDensity, EigenEnergy, Occupation, 
         EigenEnergyOnMesh, QuantumDensity)
 poisson_simulation = dd
 strain_simulation = strain
number_of_eigenstates = 12
 
Solver
{
 solver = krylovshur
}

 Physics
 {
 particle = el
 model = single_band
 }
}

This the Module declaration for electrons.

We choose a single band model with electron as particle:

Physics
 {
 particle = el
 model = single_band
 }

(equivalent to the old syntax model = conduction_band)

the number of eigenvalues to be computed is given by:

number_of_eigenstates = 12

Do not forget the links with the driftdiffusion simulation, which provides band edges:

poisson_simulation = driftdiffusion

and with elasticity simulation, which provides strain info:

strain_simulation = strain

In this way we connect Module efaschroedinger to the modules which calculate poisson and strain models. Note that this is very important, otherwise the quantum model will not have access to the information on strain and potential which are needed to write correctly the Hamiltonian of the system.

In the same way, for holes we define a quantum_hl  simulation for a k.p 6x6 calculation of valence states with:

Module efaschroedinger

{
 
 name = quantum_hl
 regions = quantum 
 plot = (ProbabilityDensity, EigenEnergy, Occupation, 
         EigenEnergyOnMesh, QuantumDensity)
 
 poisson_simulation = dd
 strain_simulation = strain
 number_of_eigenstates = 25 

 Physics
 { 
 model = 6x6
 }

}

Note that we choose a 6 band kp model with:

Physics
 {
 model = 6x6
 }

equivalent to the old syntax

model = kp
kp_model = 6x6

Run simulations

Let's now run tiberCAD to calculate strain (strain simulation), driftdiffusion (sweep for dd simulation) and quantum eigenvalues for holes and electrons (quantum_el, quantum_hl)

solve = (strain, sweep, quantum_el,quantum_hl)

tibercad InGaN_NW.tib

ATTACHMENTS

• InGaN_NW.geo
  
• InGaN_NW.msh
  
• InGaN_NW.tib