Si PN DIODE 1D

In this tutorial we will show a simple 1D Si pn diode example.
We will calculate the IV characteristic of the pn diode, by solving the Poisson equation and calculating the current with a drift-diffusion scheme.

The semi-classical transport simulation of electrons and holes is based on the driftdiffusion approximation.
Beside the electric potential, the electro-chemical potentials are used as variables such that the system of PDEs to be solved reads as follows:   where P is the electric polarization due to e.g. piezoelectric effects and R is the net recombination rate, i.e. recombination rate minus generation rate.

In order to execute correctly the example you should have the following files in the same working directory:
diode_1D.tib : input file for TiberCAD
diode_1D.msh : mesh file produced by GMSH
GMSH is an automatic three-dimensional finite element mesh generator, please refer to GMSH manual for any details.

The mesh has been obtained from the GMSH geo script file diode_1D.geo. ### MODELING OF THE DEVICE WITH GMSH

In this script, first the geometrical entities Points and Lines are defined. Then, two physical regions are defined, each associated to one of the two geometrical entities: `Physical Line("p_side")` and `Physical Line("n_side")`.

```Physical Line("p_side") = {1};
Physical Line("n _side")= {2};
```

In this way, these physical regions are made available for tiberCAD, and will be associated to a tiberCAD Region (see in the following).

Then, two Physical Points are defined, anode and cathode, and associated to the first and to the last point of our 1D device.

```Physical Point("anode") = {1};
Physical Point("cathode") = {3};
```

These points are needed to impose some boundary conditions and in this way they are made available for TiberCAD, and will be associated each to a boundary condition region (see in the following).

A non-interactive mode is also available in GMSH, without graphical user interface. For example, to mesh this first 1D tutorial in non-interactive mode, just type:

`gmsh diode.geo -1 -o diode.msh`

where diode.geo is the geometrical description of the device with GMSH syntax.

Some command line options are

• -1 , -2, -3 to perform 1D, 2D or 3D mesh generation
• -o mesh_file.msh to specify the name of the mesh file to be generated

### DEVICE STRUCTURE

Let's have a look to the input file; for further details you can refer to the reference manual.

```Device
{
meshfile = diode_1D.msh
material = Si

Region n_side
{
Doping
{
Nd = 1e18
type = donor
}
}

Region p_side
{
Doping
{
Nd = 1e18
type = acceptor
}
}
}
```

First we define the device structure: it is composed by two Regions, both constituted of the material Silicon, respectively n and p doped with a concentration of 1e18 cm-3.
With `Region n_side ` and ` Region p_side` the physical regions previously defined in GMSH are associated to the respective tiberCAD Region.

### SIMULATION MODULES

Now, we have to define the simulations which we are going to execute in our calculations.

One simulation is a particular set of equations, boundary conditions, physical parameters, solver parameters, which describes a physical problem to be solved by TiberCAD.
A valid tiberCAD simulation must belong to one of the predefined tiberCAD Modules.

TiberCAD, being a multiphysics software tool, includes several simulation models, each one describing a physical problem to be solved, e.g driftdiffusion (to solve Poisson and DriftDiffusion equations), efaschroedinger (to solve Schroedinger equation in effective mass approximation), Elasticity (to calculate macroscopical strain and ..) and others.

To create a TiberCAD simulation, we first have to declare the tiberCAD Module to which our simulation belongs:

```Module driftdiffusion
{
name = dd
```

Here, we declare that the simulation to be created will belong to the Module driftdiffusion (`Module driftdiffusion`)
Then we define the name of the TiberCAD simulation (`name = dd`) In this example, just one simulation for the driftdiffusion Module is created. In general, several tiberCAD simulations belonging to the same Module can be created; each of them must have a different name.

```Module driftdiffusion
{
name = dd
plot = (Ec, Ev, Eg, ElField, ElPotential, eQFermi, hQFermi, eDensity,
hDensity, NetRecombination, eCurrentDensity, hCurrentDensity,
CurrentDensity, ContactCurrent, eMobility, hMobility, IonizedDonors,
IonizedAcceptors)

Physics
{
particle_density
{
statistics = boltzmann
}
# use a field dependent mobility model with the
# gradient of the electrochem. potential as field
# parameter
mobility field_dependent
{
low_field_model = doping_dependent
}
# use SRH and Auger recombination
recombination (srh, auger) {}
}
}
```

Thus, in this simulation we are going to calculate Poisson and Drift-diffusion (Module driftdiffusion) for all the device.

Next, the Physics section follows.
Physics is used in general to define the physical model to be used to describe a physical property or quantity associated to the simulation.

Here, for example, we choose a Schottky-Read-Hall (SRH) and Auger models to describe recombination in our driftdiffusion calculations.

```recombination (srh, auger)
```

Also, we choose a field dependent model to describe electron and hole mobility (`keyword` `field_dependent`), with the low field mobility determined by a doping-dependent model (Masetti model, see reference manual) (`low_field_model = doping_dependent`).

```mobility field_dependent
{
low_field_model = doping_dependent
}
```
Now, we have to specify the Boundary conditions for the PDE problem represented by the simulation `dd`: they are defined in `Contact` blocks.

Here, two boundary conditions are defined, corresponding to the two contacts of our pn diode, anode and cathode.

```Contact anode
{
type = ohmic
voltage = \$Vb[0.0]
}
Contact cathode
{
type = ohmic
voltage = 0.0
}
```

Here the `anode` and `cathode` `Contacts` are associated to the corresponding `anode` and `cathode` regions defined in GMSH as Boundary condition regions (which in the present case of 1D simulation are actually just a point).

Besides, the type of Boundary condition is assigned (in this case `ohmic`) with its (initial) value, given by `voltage`.

Note that the anode voltage is expressed by the notation `\$Vb[0.0]`.
This means that the anode voltage will be given the value of the variable `Vb`, which will be specified in the `sweep` block (see below). `[0.0]` means that the default voltage value is 0.0 V.

### SWEEP BLOCK

```Module sweep
{
variable = \$Vb
solve = dd
start = 0.0
stop = 1.2
steps = 60
plot_data = true # to write data for each step
}
```

The sweep block allows us to perform the calculation of the IV characteristic of the diode: it uses the variable `\$Vb`, defined in driftdiffusion, which will assume a range of values between 0 and 1.2 V. These values will be assigned in turn to the anode contact, in order to execute a series of complete `dd` simulations for each bias condition.

### RUN SIMULATIONS

Finally, in Simulation section, we define the actual simulation to be performed, specified by

`solve = sweep`

this determines the execution of `dd` inside the `sweep` cycle;
this means that we are going to perform a calculation of Poisson and drift-diffusion for all the biases defined in the `sweep` section.

Results of the calculation will be in ascii column data format (`output_format = grace`), suited to 1D case, since it consists of ascii text files with column data.

```Simulation
{
verbose = 2
solve = sweep
resultpath = diode1D_results
output_format = grace
}
```

`tibercad diode.tib`

During the execution, the screen output shows some information about how the simulation is going:

First, about the creation of the Device regions; if there is an inconsistency between the mesh file and input file device description, here an error message is displayed and the execution terminated.

Then, the defined Modules are initialised and the selected options and available output variables are shown on screen

If all is ok, the calculation is started:
We solve: sweep
In this example only one simulation is performed, sweep, which is a set of drift-diffusion calculations.

First, the equilibrium solution is found, by solving Poisson equation:

```Ramping to sweep value `\$Vb = 0`
INFO: I will first solve a currently unsolved system (dd):

>>>>---------------------------------- (dd) ---------------
Setting max_iterations for nonlinear solver to 150
Solving equilibrium
it 1, |du| = 11.0349, |r| = 1.88523 alpha = 1
it 2, |du| = 5.30845, |r| = 0.199306 alpha = 1
it 3, |du| = 1.29314, |r| = 0.0264031 alpha = 1
it 4, |du| = 0.152896, |r| = 0.000721821 alpha = 1
it 5, |du| = 0.0034889, |r| = 5.61159e-007 alpha = 1

iterations: 6, |du| = 2.57687e-006, |r| = 3.61564e-013
Equilibrium done

......
```
Then, current continuity equations are solved together with Poisson eq., for each step of the sweep.
The voltage and current found for each contact is displayed:

```Ramping to sweep value \$Vb = 0.02
Trying \$Vb = 0.02

>>>>---------------------------------- (dd) ------------
it 1, |du| = 0.773634, |r| = 0.000488238 alpha = 1
it 2, |du| = 0.0769738, |r| = 2.33393e-005 alpha = 1
it 3, |du| = 0.00280745, |r| = 5.16137e-008 alpha = 1

iterations: 4, |du| = 4.23387e-006, |r| = 5.19017e-012

contact name: voltage: current:
anode 0.02 -2.84824e-009
cathode 0 2.69305e-009

..........................

The simulation ends correctly (hopefully...) with the Sweep value Vb = 1.2.

Ramping to sweep value \$Vb = 1.2
Trying \$Vb = 1.2

>>>>---------------------------------- (dd) --------------
it 1, |du| = 0.773634, |r| = 0.00968341 alpha = 1
it 2, |du| = 0.00292795, |r| = 2.16939e-005 alpha = 1

iterations: 3, |du| = 2.76748e-005, |r| = 7.49219e-007

contact name: voltage: current:
anode 1.2 -2526
cathode 0 2526

Solve time: 0h 0m 0.3s
<<<<------------------------------------------------------

Solve time: 0h 0m 38s
<<<<--------------------------------------------------------

Total solve time: 0h 0m 38s
Simulation finished...
Goodbye
```

### OUTPUT

Now, the output directory contain the simulation results, as defined in `output_format` and `plot`.

In dd_Vb_0_msh.dat and all the other files for each bias step, we have the output for the quantities which have been calculated, e.g. conduction and valence bands, (quasi)fermi levels, electron and hole density and mobility.

Here is, for example, the band profile obtained at the last step of calculation, for V = 2V.

The IV characteristic of the diode is contained in the file sweep_dd.dat, where the currents at both the contacts are reported.

In the file dd_Vb_0_msh.dat etc., you can find also the other quantities which have been calculated, in this case the electron, hole and total current density (Jn_x, Jp_x , J_x).