Thanks to Juan’s efforts on developing the DEVSIM. Especially after version 1.5, a lot of models can be implemented with the support of edge/element node volume based integration.

Here I wanna share an example of Avalanche Breakdown by implementing impact ionization model.

The model comes from the work in [1]. By using the constant temperature assumption, the ionization coefficient fomula can be simplified to:

`Imp_coeff = a*F_ava*exp(-(b/F_ava)^2)`

where a and b are empirical parameters and F_ava is the driving force.

The driving force can be defined as the absolute value of the gradient of quasi-fermi level:

`F_ava = |grad(EFN or EFP)|`

Since the quasi-fermi levels for electrons and holes (EFN and EFP) are already defined in devsim_bjt_example (devsim_bjt_example/simdir/physics/new_physics), we can just append their gradient and derivative definations at the end of function CreateQuasiFermiLevels:

```
edge_average_model(device=device,region=region,node_model="EFN",edge_model="dEFN",average_type='gradient')
edge_average_model(device=device,region=region,node_model="EFP",edge_model="dEFP",average_type='gradient')
for v in variables:
edge_average_model(device=device,region=region,node_model="EFN",edge_model="dEFN",derivative=v,average_type='gradient')
edge_average_model(device=device,region=region,node_model="EFP",edge_model="dEFP",derivative=v,average_type='gradient')
```

EFN:Holes and EFH:Electrons can be set to zero.

The driving forces and ionization coefficients for electrons and holes should be:

```
Fava_n = "abs(dEFN)"
Fava_p = "abs(dEFP)"
Ion_coeff_n = "a_n*Fava_n*exp(-((b_n/(Fava_n+1e-30))^2))"
Ion_coeff_p = "a_p*Fava_p*exp(-((b_p/(Fava_p+1e-30))^2))"
```

Note that both driving force and current are along the edge. We assume they were constant in the area corresponding to the “EdgeNodeVolume” defined in [2]. Thus the impact ionization rate in this area can be calculated as:

`Imp_rate = "(Ion_coeff_n*(abs(Jn))+Ion_coeff_p*(abs(Jp)))/q"`

which is also an edge model. The above physical quantities and their derivatives to the variables should be defined before current continuity equation:

```
CreateEdgeModel(device,region,"Fava_n",Fava_n)
CreateEdgeModel(device,region,"Fava_p",Fava_p)
CreateEdgeModel(device,region,"Ion_coeff_n",Ion_coeff_n)
CreateEdgeModel(device,region,"Ion_coeff_p",Ion_coeff_p)
CreateEdgeModel(device,region,"qGn_imp","q * %s" % Imp_rate)
CreateEdgeModel(device,region,"qGp_imp","-q * %s" % Imp_rate)
for v in variables:
CreateEdgeModelDerivatives(device,region,"Fava_n",Fava_n,v)
CreateEdgeModelDerivatives(device,region,"Fava_p",Fava_p,v)
CreateEdgeModelDerivatives(device,region,"Ion_coeff_n",Ion_coeff_n,v)
CreateEdgeModelDerivatives(device,region,"Ion_coeff_p",Ion_coeff_p,v)
CreateEdgeModelDerivatives(device,region,"qGn_imp","q * %s" % Imp_rate,v)
edge_model(device=device,region=region,name="qGp_imp:%s" % v,equation="-qGn_imp:%s" % v)
```

Next, put the “qGn_imp” and “qGp_imp” to field *edge_volume_model* in ECE and HCE, respectively.

Then you can insert a big resistor in series with your device and sweep the votlage to simulate the avalanche breakdown.

Ref:

[1] Okuto, Y., and C. R. Crowell. “Threshold energy effect on avalanche breakdown voltage in semiconductor junctions.” *Solid-State Electronics* 18.2 (1975): 161-168.

[2] Sanchez, Juan, and Qiusong Chen. “Element Edge Based Discretization for TCAD Device Simulation.” (2021).