11institutetext: Physics Division, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA22institutetext: Department of Physics and Astronomy,
Texas Tech University, Lubbock, TX 79409, USA33institutetext: Electronic
Engineering Dept., School of Engineering, University of Sevilla, 41092
Spain44institutetext: Nikhef National Institute for Subatomic Physics, Science
Park 105, 1098 XG Amsterdam, Netherlands55institutetext: Fermilab, P.O. Box
500, Batavia Il USA66institutetext: INFN Sezione di Perugia, Perugia,
Italy77institutetext: Consiglio Nazionale delle Ricerche - Istituto Officina
dei Materiali, Perugia, Italy88institutetext: Université de Paris &
Laboratoire Astroparticule et Cosmologie, Paris, France99institutetext:
Department of Physics and Astronomy, University of New Mexico, Albuquerque, NM
87131, USA1111institutetext: Engineering Department, University of Perugia,
Perugia, Italy1212institutetext: Institute of High Energy Physics, Chinese
Academy of Sciences, Beijing 100049, China

# Simulations of Silicon Radiation Detectors for
High Energy Physics Experiments

B. Nachman (ed.), 2 T. Peltola (ed.) 6, 7 P. Asenov 8 M. Bomben 5 R. Lipton 6,
7 F. Moscatelli 9 E. A. Narayanan 3 F. R. Palomo 6, 11 D. Passeri 9 S. Seidel
12 X. Shi 4 J. Sonneveld<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Silicon radiation detectors are an integral component of current and planned
collider experiments in high energy physics. Simulations of these detectors
are essential for deciding operational configurations, for performing precise
data analysis, and for developing future detectors. In this white paper, we
briefly review the existing tools and discuss challenges for the future that
will require research and development to be able to cope with the foreseen
extreme radiation environments of the High Luminosity runs of the Large Hadron
Collider and future hadron colliders like FCC-hh and SPPC.

Submitted to the Proceedings of the US Community Study

on the Future of Particle Physics (Snowmass 2021)

## Executive Summary

There are currently a variety of tools available for simulating the properties
of silicon sensors before and after irradiation. These tools include finite
element methods for device properties, dedicated annealing models, and
testbeam/full detector system models. There is no single model that can
describe all of the necessary physics. Most of these models are either fully
or partially developed by HEP scientists and while there are many open-source
tools, the most precise device property simulations rely on expensive,
proprietary software.

The development of these simulations happens inside existing experimental
collaborations and within the RD50 Collaboration at CERN. RD50 is essential
for model research and development and provides an important forum for inter-
collaboration exchange.

While existing approaches are able to describe many aspects of signal
formation in silicon devices, even after irradiation and annealing, there is
significant research and development (R&D) required to improve the accuracy
and precision of these models and to be able to handle new devices (e.g. for
timing) and the extreme fluences of future colliders. The US particle physics
community can play a key role in this R&D program, but it will require
resources for training, software, testbeam, and personnel.

For example, there is a great need for (1) a unified microscopic model of
sensor charge collection, radiation damage, and annealing (no model can
currently do all three), (2) radiation damage models (for leakage current,
depletion voltage, charge collection) with uncertainties (and a database of
such models), and (3) a measurement program to determine damage factors and
uncertainties for particle types and energies relevant for current and future
colliders.

##  1 Introduction

Silicon (Si) devices are now an integral component of collider-based high
energy physics (HEP) experiments. These detectors are being used as tracking
detectors with precise space (and now time) measurements, for planned
calorimeter detectors, and for beam monitoring. Detailed Si simulations are
essential tools for running experiments and for developing experiments of the
future. Simulations are used to model existing detectors in order to determine
operational parameters like temperature and high voltage. They are also needed
to accurately describe the experiment for data analysis. Simulations
additionally are also used for device structure optimization and other
detector design tasks. Silicon simulation at the Large Hadron Collider (LHC)
and future colliders requires dedicated efforts because of the unique devices
needed by our community and because of the extreme radiation environments to
which these devices are exposed. A detailed review of radiation effects in the
context of LHC experiments can be found at Ref. Dawson:2764325 .

Multiple types of Si simulations are required to serve all of the needs of
detector operation, analysis, and design. At the lowest level, device
simulations are used to predict the local electric field and signal formation.
Often these tools do not model the time dependence of the detector response,
which is influenced by thermal annealing. Testbeam and full detector systems
also must be modeled in order to compare with data from colliders.

At a microscopic level, the accumulation of radiation-induced defects
progressively degrades various operational properties of a silicon sensor both
inside the Si-bulk and at the surface on the Si/SiO2-interface. These
properties are presently addressed by several defect models that apply 2–5
acceptor and donor type trap levels in the Si-substrate, and fixed oxide
charge density ($N_{\textrm{f}}$) with/without 1–2 acceptor and donor type
interface-trap ($N_{\textrm{it}}$) levels at the surface. At a macroscopic
level, silicon devices are characterized by their signal formation, leakage
current, depletion voltage, etc. The signal formation is then responsible for
higher-level quantities like track reconstruction and physics analysis with
charged-particle tracks.

The goal of this report is to consider Si simulations in the context of the US
Community Study on the Future of Particle Physics (Snowmass 2021). The
following discussion is divided into two parts, introduction of existing
simulation tools and radiation defect models and then the challenges faced in
their further development for future requirements. In the first part, after
describing the mathematical models applied in the numerical Technology
Computer-Aided Design (TCAD)-simulations in Section A, the currently utilized
tools for the per-sensor simulations, as well as for the testbeam-telescope-
setup and full-detector-system simulations are introduced in Section 2. Also
presented in Section 2, and based largely on the results summarized in the
Chapter 7 of the recent CERN Yellow Report Dawson:2764325 , are the various
radiation defect models for bulk and surface damage or their combination, that
have been developed in the experiments of the LHC. Additionally, speedup
investigations of large-mesh-size ($\sim 10^{5}-10^{6}$ mesh-points) TCAD
devices and TCAD’s process simulation option are discussed in Sections 2.1.2
and 2.2, respectively. The second part of the paper in Section 3 addresses
both the observed and anticipated requirements and challenges on the
simulation frameworks and the further development of radiation defect models
toward more unified approaches that are able to account for the extreme
radiation environments in the future colliders. Potential issues with
proprietary software and computational demands of large mesh-size devices are
discussed in Sections 3.1 and 3.2, respectively, while the development needs
of simulation chains between packages and the sources of uncertainties to be
taken into account in the TCAD modeling are addressed in Sections 3.3 and 3.4,
respectively. Finally, the approaches to upgrading and unifying existing
defect models in the radiation fluence region beyond $1\times 10^{16}$
n${}_{\textrm{\small eq}}/$cm2, expected for future colliders, are considered
in Sections 3.5, 3.6 and 3.8.

##  2 Existing Tools

###  2.1 TCAD Simulations for Detector Properties

Radiation levels above $\sim
10^{13}~{}\textrm{n}_{\textrm{eq}}/\textrm{cm}^{2}$ introduce observable
damage to the crystal structure of a silicon sensor. Fluences beyond $1\times
10^{14}~{}\textrm{n}_{\textrm{eq}}/\textrm{cm}^{2}$ lead eventually to a
significant degradation of the detector performance. In the immense radiation
environment of the HL-LHC, defects are introduced both in the silicon
substrate (bulk or displacement damage) and in the SiO2 passivation layer,
that affect the sensor performance via the interface with the silicon (surface
damage). The multitude of observed defect levels after irradiations with
hadrons or higher energy leptons Radu2015 , set up a broad parameter space
that is not practical to model and tune. Thus, minimized sets of defects
constituting various effective defect models have been applied as the approach
for the simulations of irradiated silicon detectors. Most of these defect
models have been developed for the commercial finite-element TCAD
software111See Appendix A for the governing equations. frameworks Synopsys
Sentaurus™222http://www.synopsys.com and Silvaco
Atlas™333http://www.silvaco.com.

The modeling of the main macroscopic effects of the bulk damage on high-
resistivity Si-sensors irradiated by hadrons–the change of the effective
doping concentration $N_{\textrm{eff}}$ (resulting in modified electric field
distribution and change of operational voltage), the increase in the leakage
current proportional to the fluence and the degradation of charge collection
efficiency (CCE)–includes two-level Eremin2002 ; Chiochia2005 ; Eremin2011 ;
Verbit2012 ; Eber2013 , three-level Petasecca2006 ; Pennicard2008 ;
Passeri2016 ; FOLKESTAD201794 , as well as a five-level models schwandt2019new
. Figure 1 demonstrates a comparison between measured transient signals
(Figure 1a) and CCE (Figure 1b) of neutron irradiated sensors and simulated
results produced by a two-level defect model for neutrons.

Simulating the effects of the surface damage on Si/SiO2-interface irradiated
by charged hadrons, gammas, or X-rays –like the change of inter-electrode
resistance ($R_{\textrm{int}}$) and capacitance ($C_{\textrm{int}}$), modified
electric fields at implant edges and charge-injection position dependence of
CCE (CCE($x$))– involves approaches, where surface damage is approximated
solely by the fixed oxide charge density ($N_{\textrm{f}}$) Verzellesi2000 ;
Piemonte2006 ; Unno2011  or by three interface traps ($N_{\textrm{it}}$) with
parameters matching the measurements of X-ray irradiated MOS-capacitors
Zhang2012 . Results from a further model that combines the two approaches is
presented in Figure 2, where the characteristics of the $CV$-curves of a
$\gamma$-irradiated MOS-capacitor are only reproduced by simulation by
including in addition to $N_{\textrm{f}}$ both acceptor and donor-type deep
$N_{\textrm{it}}$ ($N_{\textrm{it,acc,don}}$) at the Si/SiO2-interface (Figure
2a). By using both $N_{\textrm{f}}$ and $N_{\textrm{it,acc,don}}$ as an input
for an $R_{\textrm{int}}$-simulation with typical CMS High Granularity
Calorimeter (HGCAL) isolation implant parameters in Figure 2b, the
experimentally observed $R_{\textrm{int}}$-values are quantitatively
reproduced (atoll and common p-stops with STD $N_{\textrm{ps}}$). By excluding
$N_{\textrm{it,acc,don}}$ from the simulation, the pads become either shorted
(STD $N_{\textrm{ps}}$, $N_{\textrm{it}}=0$) for all voltages or reach
isolation only above 500 V of reverse bias voltage for an extreme value of
p-stop peak doping ($5\times\textrm{STD}$ $N_{\textrm{ps}}$,
$N_{\textrm{it}}=0$).

First steps towards a unified bulk and surface defect model include an
approach where two bulk defect levels are augmented by an acceptor-type
$N_{\textrm{it}}$ with 2 $\upmu$m depth distribution from the surface
(Sentaurus Peltola2014 ; Peltola2015 ), and another approach where two bulk
defect levels are complemented by two acceptor-type $N_{\textrm{it}}$ (Silvaco
Dalal2014b ). More details on the development of the defect models are
described in e.g., Peltola2015r .

(a)

(b)

Figure 1: Examples of bulk-damage modeling Peltola2020 . (a) Measured and
simulated IR-laser induced transient signals for 300-$\upmu$m-thick $p$-on-$n$
pad-sensor (300N) after neutron irradiation to the fluence of $(6.1\pm
0.5)\times 10^{14}~{}\textrm{n}_{\textrm{eq}}/\textrm{cm}^{2}$. The simulation
applied neutron defect model Eber2013  with a fluence of $6.0\times
10^{14}~{}\textrm{n}_{\textrm{eq}}/\textrm{cm}^{2}$. The sensor parameters
used in the simulation were extracted from $CV/IV$-measurements. (b) Measured
and simulated evolution of charge collection efficiency (CCE) with voltage and
1-MeV equivalent neutron fluence for 300-$\upmu\textrm{m}$ thick sensors at
-30 ∘C. Simulated results use dotted, dashed and solid curves for 600 V, 800 V
and 1 kV, respectively.

(a)

(b)

Figure 2: Examples of surface-damage modeling. (a) Measured and simulated
$CV$-curves at $T=\textrm{RT}$ and $f=2$ kHz of a $\gamma$-irradiated
200-$\upmu$m-thick $n$-bulk MOS-capacitor with $t_{\textrm{ox}}=700$ nm for
the dose of $7.0\pm 0.4$ kGy. Measurements included voltages starting from
both inversion (Inv.) and accumulation (Acc.) regions. (b) Inter-pad
resistivity at -20∘ C for individual (atoll) and common p-stop isolation
implants with HGCAL-parameters, simulated by applying $N_{\textrm{f}}$ and
$N_{\textrm{it}}$-parameters tuned from MOS-capacitor measurements and
simulations for the dose of 23.5 kGy. The 100 M$\Omega\cdot$cm (black dashed
line) represents conservative estimation for high level of pad isolation. STD
$N_{\textrm{ps}}$ = standard value of p-stop peak doping.

Another set of approaches to unify the surface and bulk defect models with
common characteristics is known as the Perugia model (named after the city of
the group which developed it). It takes into account the fact that in
principle, the radiation damage in silicon sensors is caused by two different
factors: the ionizing energy loss and the non-ionizing energy loss. Although
the ionizing energy loss is important for the signal formation and is usually
reversible, it can cause irreversible effects on the oxide by introducing
positive oxide charge in the SiO2, by increasing the number of bulk oxide
traps and by increasing the number of interface traps. These effects are
referred to as surface damage in the model, and they influence the operation
of segmented silicon sensors with respect to the inter-electrode isolation,
the breakdown voltage and the charge collection efficiency in real-life
experiments. The non-ionizing energy loss is responsible for the introduction
of defects into the silicon lattice through the displacement of crystal atoms,
usually produced due to the impact of high-momentum particles. These impacts
lead to point and cluster defect generation and hence to the introduction of
deep-level trap states which act like generation-recombination centers. Non-
ionizing effects are referred to as bulk damage in the Perugia model and on a
macroscopic scale they are responsible for the increase of the leakage current
in silicon sensors, the changes in the effective space charge concentration
and the charge collection efficiency Moscatelli:2017sps .

In the original Perugia model Passeri2016 ; Moscatelli:2017sps , a three level
model is presented for simulating the bulk damage effects for $n$-type and
$p$-type substrates. In order to incorporate also the surface damage effects
an extension of the model has been made referred to as the Perugia surface
model 2019 Morozzi:2020yxm ; Morozzi:2021kwo . The surface damage effects can
be mainly described by three parameters: the oxide charge density
($N_{\textrm{ox}}$) and the two interface trap states ($N_{\textrm{it}}$) for
donors and acceptors. The values of the above quantities can be extracted from
high-frequency and quasi-static $CV$-measurements and then used as inputs in
the simulation. Macroscopic factors, such as the capacitance or the leakage
current, are calculated as outputs. This version is able to reproduce the
radiation damage macroscopic effects up to the order of
$10^{15}\,$n${}_{\textrm{eq}}$/cm${}^{2}\,$ 1-MeV neutron equivalent fluences
and X-ray photon doses of up to 10 Mrad. The Perugia surface model 2019 damage
modeling scheme, is fully implemented within the Sentaurus TCAD environment.
Its main advantages are the dependence on a limited number of parameter
relevant for physics and the versatility of the simulation approach. Up to
now, it has already been compared with experimental measurements performed on
Gate-Controlled Diodes (GCD) (Figure 3) and Metal-Oxide-Semiconductor (MOS)
capacitors at Fondazione Bruno Kessler (FBK) in Italy, Hamamatsu Photonics
K.K. (HPK) in Japan and Infineon Technologies (IFX) in Austria on both
$n$-type and $p$-type substrates. Furthermore, an updated radiation damage
model (called New University of Perugia TCAD model) has been fully implemented
within the simulation environment, to have a predictive insight into the
electrical behavior and the charge collection properties of Low Gain Avalanche
Detectors (LGAD) detectors. By coupling the New University of Perugia TCAD
model with an analytical model that describes the mechanism of acceptor
removal in the multiplication layer of the LGAD, it has been possible to
reproduce experimental data with high accuracy (Figure 4) Croci:2022ygd .
Thanks to the intrinsic multiplication of the charge within the LGADs, it is
possible to improve the signal-to-noise ratio (SNR), thus limiting its drastic
reduction with fluence, and hence a special attention has been given to the
choice of the avalanche model, which allows the simulation findings to better
fit with experimental data. The Massey avalanche model was chosen as the best
option since it is the only one that does not underestimate the breakdown
voltage of the device and leads to successful comparison with the experimental
data Massey .

This last version of the Perugia model is able to reproduce the behavior of
3D-devices as well.

Figure 3: $I$-$V$ characteristics of $n$-on-$p$ HPK GCD: comparison between
the measured and simulated diode current as a function of gate voltage (with
$p$-spray isolation layer). Measurements in solid line and simulations in
dashed line. The model used in the simulation is Perugia surface model 2019.

(a)

(b)

Figure 4: (a) $I$-$V$ and (b) $C$-$V$ simulated curves compared with
experimental data before and after irradiation for a Low-Gain Avalanche Diode
(LGAD) produced by FBK. The sensor area is $1\,$mm2, the thickness is
55$\,\upmu$m and the temperature is 300 K. In the simulation, the New
University of Perugia TCAD model was coupled with an analytic model developed
by the Torino group that describes the mechanism of acceptor removal in the
multiplication layer.

####  2.1.1 TCAD Alternatives

There are simpler software alternatives to TCAD like TRACS, tracs ,
Weightfield2, weightfield  or KDetSim444http://kdetsim.org kdetsim1 . All of
them solve the Poisson equation (7) for a charge distribution
$N_{\textrm{eff}}$. Instead of calculating $N_{\textrm{eff}}$ from the full
set of continuity equations, carrier drift is done in a static electric field
and calculated by steps and the detector signal (induced current) comes from
the Shockley-Ramo theorem, shockley ; ramo :

$I_{i}=-q\overrightarrow{v}(x)\cdot\overrightarrow{E}_{w}(x)\,,$ (1)

where $x$ is the carrier position, $v(x)$ is its drift velocity that depends
on the applied electric field $E(x)$ and $E_{w}(x)$ is the weighting field for
the particular $i$-electrode. This lightweight approach does not substitute a
TCAD software but it is appropriate, for example, to design a full TCAD
simulation and for obtaining approximate results much faster than with any
TCAD suite, which would be more precise and capable but also more demanding in
computational resources. In any case, none of the lightweight alternatives are
currently suited to calculate trap dynamics.

####  2.1.2 Simulations of Large-Mesh Devices

In the Synopsys Sentaurus TCAD software framework, the device structures can
be generated in both 2D or 3D. Sensors that have negligible contribution to
the weighting field from the third dimension, i.e. diodes and strip sensors,
can be accurately modeled in 2D, and extended to the dimensions of a real
device by an appropriate factor. The 3D device structures, although requiring
much computing time and processing capacity, are mandatory for reliable
simulations of the pixel and 3D-columnar sensors, presented in Figure 5. In a
typical parameter scan the number of nodes in the simulation can be 10–20,
while the simulation of a single node for a 3D-structure with large enough
mesh-size to reproduce a realistic device (e.g. smooth implant shapes at the
electrodes) can be in the order of days in a standard multi-core processor.
Therefore, considering different options for the optimization of the
simulation times is highly motivated.

Speedup investigation of a 3D-device with about $10^{6}$ mesh-points in Figure
6 shows, that shared-memory parallelization benefits simulation execution time
by a factor of five when the number of threads is increased from one to match
the number of physical cores in the processor. A further factor of about three
is gained when the computation method of the physics equations is changed from
a direct linear solver (optimal for 2D-simulations) to an iterative linear
solver (ILS) with tuned parameters to avoid compromising the accuracy of the
simulation results. Weakly supported distributed processing (cluster
computing) in Synopsys TCAD is possible to be compensated in a cluster by
generating copies of the simulation project that are dedicated for each node
of the parameter space and running these simultaneously.

(a)

(b)

(c)

Figure 5: 3D-sensor structures designed in Sentaurus TCAD. (a) Electric fields
in a $50\times 50$ $\mu\textrm{m}^{2}$ $n$-on-$p$ pixel sensor for CMS
Tracker. (b) Electric fields in a $p$-on-$p$ MediPix sensor Peltola2016p  with
a 55 $\mu\textrm{m}$ pitch. (c) Sliced view of a double column-double sided 3D
$n$-on-$p$ sensor with a p-stop isolation. Aluminum and oxide have been
stripped from the surface. Figure 6: Speedup by shared-memory parallelization
in a 3D-detector structure with about $10^{6}$ mesh-points, where the physics
equations are solved. Theoretical speedup courtesy of Synopsys support.

###  2.2 Process Simulations for Device Development

Commercial TCAD packages are capable of full simulations of device
fabrication, including epitaxy, implantation, annealing, deposition and
oxidation. The accuracy and detail provided by these simulations can be
invaluable in the development of new sensor technologies or in understanding
the behaviour of existing devices. For devices such as LGADs, an accurate
model of the implant doping density and profile is crucial to understanding
and predicting performance. This is possible with process simulation as the
simulation includes detailed models of the implantation and annealing
including crystal channeling, implantation damage, and dopant activation. The
process model can also provide a link between commercial foundries and
instrumentation developers. Figure 7 shows the result of the process
simulation for 8-inch AC-coupled wafer development. The wafers were fabricated
by NHanced/Novati with process flow input from LBNL and SLAC and the Silvaco
process simulation from Fermilab Alyari:2020pet .

Figure 7: Results of the process simulation of the strip end region of an AC
coupled detector including implants, AC coupling layers, and polysilicon
resistors.

###  2.3 Hamburg Model Simulation of Silicon Radiation Damage

####  2.3.1 Modeling Silicon Leakage Current: the Hamburg Model

One of the best-characterized methods for monitoring silicon radiation damage
is based on measuring the sensor leakage current. The most widely used-model
for predicting the leakage current is the Hamburg Model moll-thesis . A linear
relationship between particle fluence and leakage current applies to silicon
sensors. The effective fluence, $\Phi_{\mathrm{eq}}$, is the number of $1$-MeV
neutrons applied to a sensor of surface area $1$ cm2 that cause damage
equivalent to that of all particles that went through the sensor. The linear
relationship between effective fluence and $\Delta I$, the change in leakage
current at fluence $\Phi_{\mathrm{eq}}$ relative to the value before
irradiation, is given as

$\Delta I=\alpha\cdot\Phi_{\mathrm{eq}}\cdot V\,,$ (2)

where $V$ is the depleted volume of the silicon sensor and $\alpha$ is the
current-related damage coefficient Lindstrom:2002gb  which is a function of
annealing time and temperature. The observed leakage current and predictions
are typically normalized to the same temperature for direct comparison. The
normalization is made using the known leakage current dependence on
temperature and the effective silicon band gap energy, $E_{\mathrm{eff}}$
Chilingarov_2013 . The change in the leakage current is predicted using the
Hamburg Model by inputting the complex radiation fields in the detector in
fluence simulated by packages such as Fluka Baranov:2005ewa ;
Battistoni:2007zzb  and Geant4 GEANT4:2002zbu ; Allison:2006ve ;
Allison:2016lfl , as given by the following equation:

$\Delta I=\left(\Phi_{\mathrm{eq}}/L_{\mathrm{int}}\right)\cdot
V\cdot\sum_{i=1}^{n}L_{\mathrm{int},i}\cdot\left[\alpha_{I}\exp\left(-\sum_{j=i}^{n}\frac{t_{j}}{\tau(T_{j})}\right)+\alpha^{*}_{0}-\beta\log\left(\sum_{j=i}^{n}\frac{\Theta(T_{j})\cdot
t_{j}}{t_{0}}\right)\right],$ (3)

where $L_{\mathrm{int}}$ is the integrated luminosity, $t_{i}$ is the time,
and $T_{i}$ is the temperature in period $i$. The first sum is over all time
periods and the two sums inside the exponential and logarithmic functions are
over the time between the irradiation in time period $i$ and the time of the
measurement. The other symbols in Eq. 3 are $t_{0}=1$ min, $V=$ depleted
sensor volume (in $\mathrm{cm}^{3}$), $\alpha_{I}=(1.23\pm 0.06)\times
10^{-17}$ A/cm, $\tau$ follows an Arrhenius equation
$\tau^{-1}=(1.2^{+5.3}_{-1.0})\times
10^{13}~{}\mathrm{s}^{-1}\times\exp^{(-1.11\pm 0.05)/k_{\mathrm{B}}T}$ (where
the units of $k_{\mathrm{B}}T$ are eV), $\alpha^{*}_{0}=7.07\times
10^{-17}$A/cm, and $\beta=(3.29\pm 0.18)\times 10^{-18}$ A/cm. The temperature
scaling function $\Theta(T)$ is defined by

$\Theta(T)=\exp\left[\frac{-E_{\mathrm{eff}}}{k_{\mathrm{B}}}\left(\frac{1}{T}-\frac{1}{T_{\mathrm{R}}}\right)\right],$
(4)

where $E_{\mathrm{eff}}=1.21$ eV is used for the effective silicon band gap
energy, $k_{\mathrm{B}}$ is the Boltzmann constant, and
$T_{\mathrm{R}}=21^{\mathrm{o}}$C. Once the simulation is complete, a
temperature correction to the full simulation is made so that it may be
compared to the leakage current data.

The leakage current predictions ATLAS:2021gld  for ATLAS outer barrel layers -
B-Layer, Layer-1, and Layer-2 are shown in Figure 8. The average of measured
leakage currents of all layers are shown for the comparison. It is clear from
the plot that the Hamburg Model correctly describes the leakage current
evolution of the ATLAS Pixel detector, after scaling to the data.

Figure 8: Average measured leakage current of a representative sample of
modules on the B-Layer, Layer-1 and Layer-2 over the full period of operation.
The scaled prediction from the Hamburg Model is also shown. The bands include
uncertainties on the measurement ATLAS:2021gld .

####  2.3.2 Modeling Depletion Voltage and Charge Collection Efficiency

The full depletion voltage $V_{\mathrm{dep}}$ of a sensor depends linearly on
the doping concentration $N_{\mathrm{eff}}$ and the square of the thickness
$d$ of the detector (wafer thickness for a planar sensor or electrode-to-
electrode distance for a 3D sensor), and can be written as

$V_{\mathrm{{dep}}}=\frac{\mid N_{\mathrm{eff}}\mid
d^{2}q}{2\epsilon\epsilon_{0}},$ (5)

where $q$ is the charge of the electron, $\epsilon$ is the dielectric
constant, and $\epsilon_{0}$ is the vacuum permittivity. The doping
concentration is affected by the fluence, and the change in doping
concentration $\Delta N_{\mathrm{eff}}$ can be expressed as

$\Delta N_{\mathrm{eff}}=N_{\mathrm{A}}+N_{\mathrm{Y}}+N_{\mathrm{C}},$ (6)

where $N_{\mathrm{A}}$, $N_{\mathrm{Y}}$, and $N_{\mathrm{C}}$ are the changes
in the doping concentration due to beneficial annealing, reverse annealing,
and stable damage respectively. The Hamburg Model555This model is not directly
related to the Hamburg Model for leakage current, even though they share the
same name. predicts the depletion voltage of a sensor when it is exposed to
radiation causing non-ionizing energy loss (NIEL) for a certain duration, by
using the change in the doping concentration described by Eq. 6. However, the
prediction requires a set of parameters called introduction rates which depend
on the material of the detector and the type of radiation the detector
receives. The existing predictions of depletion voltage use introduction rates
calculated using data recorded at relatively lower fluences moll-thesis  and
are not suitable for predictions of depletion voltage of detectors designed
for future conditions. As an illustration of this limitation, a comparison
Dawson:2764325  of depletion voltage prediction to measured depletion voltage
of the ATLAS Insertable B-Layer is shown in Figure 9.

Figure 9: Predicted Dawson:2764325  depletion voltage of the IBL according to
the Hamburg model as a function of time from installation until the end of Run
2. Circular and square points are measured depletion voltages.

As of now, the full ATLAS Pixel sensor depletion voltages recorded in Run 2
are accessible, so the Hamburg Model can be used in a fit to the data to
extract these parameters; this information can be projected to higher
fluences, and used for developments of future particle physics detectors.

An important consideration is the loss of charge collection efficiency (CCE).
CCE loss due to radiation becomes the leading failure mode at future
colliders. In the current implementation of the Hamburg Model, prediction of
charge collection is not included. Predicting the charge collection efficiency
loss of subdetectors, based on received fluence and annealing time, is
essential to ensure an efficiently functioning detector.

###  2.4 RASER – SiC detector simulation

To explore the potential capability of wide band gap semiconductor especially
silicon carbide (SiC) to be used in high energy physics experiments, a fast
simulation software – RASER (RAdiation
SEmiconductoR)666https://pypi.org/project/raser has been developed as a python
module. In RASER, the electric field and weighting potential are calculated by
FEniCS FEniCS  – an open-source computing platform for solving partial
differential equations. The simulation of electrons from 90Sr and the
deposited energy in SiC are carried out by Geant4\. The current induced by
electron–hole ($e$–$h$) pairs moving is calculated by Shockley-Ramo’s theorem
shockley ; ramo . The readout electronics used a simplified current amplifier
as in weightfield2 .

RASER has been validated by comparing measurement and simulation of time-
resolution of planar 4H-SiC detector designed by Nanjing University (NJU) in
China. The cross section of the SiC detector is shown in Figure 10. The size
of the detector is 5 mm$~{}\times$ 5 mm, and the upper and lower electrodes
are ohmic contacts. The detector has a 100 µm high resistance active 4H-SiC
epitaxial layer and a 350 µm substrate. The comparison between wave-forms
obtained from measurement, RASER, and TCAD simulation shows good consistency
as in Figure 11. Using the same parameters, the time resolution of measurement
and simulation are (94 $\pm$ 1) ps and (73 $\pm$ 1) ps, respectively sicpin .

Figure 10: Cross section of $\rm 5~{}mm\times 5~{}mm$ planar 4H-SiC detector
designed by NJU. It has $\rm 100~{}\upmu m$ 4H-SiC active epitaxy layer and a
$\rm 350~{}\upmu m$ substrate. Figure from sicpin .  Figure 11: Waveform
comparison: measured data (red) and RASER simulation (green). Figure from
sicpin .

RASER has also been used to predict the timing performance of 3D-SiC with a
potential of more radiation hard that could operate at room temperature. The
first structure of 3D 4H-SiC detector prototype studied is shown in Figure
12a, where the $p^{+}$ electrode column that penetrates the entire high-
resistance substrate perpendicular to the detector surface is surrounded by
six $n^{+}$ electrode columns. The simulated size of the 3D 4H-SiC detector is
1 cm $\times$ 1 cm, and the thickness is 350 µm, and the effective
concentration of n-type substrate is set to 1 $\times$ 1013 cm$-$3. The radius
of electrode column is 50 µm, and the column spacing, which is the distance
between the center of $n^{+}$ and $p^{+}$ electrode columns is 150 µm. The
simulated electric field is shown in Figure 12b. The 3D-SiC detector is placed
with a 34 ps LGAD under the 90Sr source where the emitted two electrons are
simulated by Geant4 as shown in Figure 12c.

The simulated time resolution reached around 35 ps with bias voltage of 300 V
(Figure 13a) while remaining constant with temperature even with values up to
+127∘ C (Figure 13b). Simulation of 5 electrodes configuration (3D-4H-SiC-5E)
reveals up to 25 ps time resolution can be reached 3D-SiC .

(a)

(b)

(c)

Figure 12:  (a) 3D 4H-SiC detector schematic diagram. (b) The electric field
distribution of x–y cross section for the 3D 4H-SiC detector with 500 V bias
voltage simulated with RASER. (c) Geant4 simulated two electrons passing
through aluminum foil, the LGAD, a printed circuit board (PCB), and the 3D
4H-SiC detector. The tracks of electrons and the deposited energy with a step
less than 1 µm were recorded. Figure from 3D-SiC .

(a)

(b)

Figure 13:  (a) Time resolution versus voltage. (b) Time resolution versus
temperature. Figure from 3D-SiC .

Future development of RASER will include the LGAD-SiC as well as the pixel-
timing performance for SiC devices.

###  2.5 Additional Capabilities

In addition to core semiconductor modelling commercial packages often have
additional functionality for specialized applications. These include
optoelectronic simulation, lattice heating simulators, single event effects,
VCSEL and LED simulation, quantum effects simulation and a number of other
options. These are typically individually licensed.

An example of such a simulation in HEP is single event burnout. Test beam
studies of irradiated LGADs showed a high rate of beam-induced breakdown
resulting in physical destruction of the device. The Silvaco Giga lattice
heating toll was used in combination with avalanche models to study the
systematics of these effects, including linear energy transfer threshold and
the effect of the field strength on breakdown margin 2010270 . Figure 14 shows
transient curves of temperature and cathode current for a 50-$\upmu$m-thick
diode operated at 625 V struck by particles with ionization between 32 and 140
MIPs.

(a)

(b)

Figure 14: Transient simulation of junction temperature (a) and cathode
current (b) for a 50-$\upmu$m-thick diode operated at 625 V exposed to heavily
ionizing particles between 32 and 140 MIPs. The threshold for burnout is
between 64 and 100 MIPs. There is about a 100 microsecond delay between impact
and the destructive heating of the sensor.

###  2.6 Testbeam

Sensor simulations from commercial packages like TCAD are less often used for
Monte Carlo simulations of drift and diffusion of carriers or, for example,
simulations of beam tests at facilities with charged particle beams. For such
purposes, specialized Monte Carlo software is used. These Monte Carlo
simulation packages can use sensor simulations results such as the sensor
electric field as input; see for example allpix2withTCAD .

One example of a software package for drift and diffusion modeling of carriers
in silicon sensors is PixelAV Chiochia2005 ; Swartz2006212 ; Swartz:687440 ,
originally created for interpretation of data taken with both unirradiated and
irradiated sensors at pion test beams. This package has accurate models of
charge deposition and transportation and can simulate charge drift under
magnetic fields. It also models radiation damage with the trapping of charges.
It uses input electric fields from the TCAD simulation software. PixelAV is
now used as standalone software in the CMS template-based hit reconstruction
software that includes radiation damage simulation Swartz:2007zz ;
Chiochia:2008yga ; CMS:2014pgm . In the templates created with PixelAV provide
both corrections to the hit position and the cluster charge.

Packages similar to PixelAV are Allpix allpix  and KDetSim. KDetSim is based
on the CERN ROOT software package. It is fast, allows for simulation of larger
volumes, and allows for iterative approaches. All packages are free and open
source.

Another free and open source package that includes experimental setups such as
beam telescopes and their material effects like multiple scattering and
nuclear interactions based on Geant4 is Allpix2 allpix2 . This includes energy
deposition based on Geant4, drift and diffusion for charge propagation, and
models charge digitization. A visualization of a setup provided by the
software package is shown in Fig. 15.

Figure 15: Visualization of a simulated detector telescope setup with 7
planes. in Allpix2. The beam is incident from the right, and the colored lines
are the primary and secondary particles propagated through the setup by the
simulation. Figure from allpix2 .

The Allpix2 framework can also provide observables for comparison with test
beam data such as particle clusters, tracks, residuals, and cluster charge. An
example is shown in Fig. 16.

(a)

(b)

Figure 16: (a) Residuals and (b) cluster charge distribution simulated with
the Allpix2 software package compared to data. The cluster charge is for a 50
$\upmu$m device under test (DUT). The residuals are the difference between the
reconstructed telescope track position in the DUT and the true hit position of
the simulated particle. Figures from allpix2 .

The packages now also include different charge carrier mobility models and the
possibility to load doping profiles for doping-dependent charge carrier
lifetime calculations spannagel2021allpix2 , and foresees to include trapping
of charges.

###  2.7 Full detector systems

Within large HEP collaborations, such as those of the LHC-experiments, silicon
detector simulation is typically performed using dedicated, proprietary777This
is changing - many collaborations now make their source code public. software
written within the experiment’s software framework (for example
Athena888https://gitlab.cern.ch/atlas/athena,
Gauss999https://gitlab.cern.ch/lhcb/Gauss, and
CMSSW101010https://github.com/cms-sw/cmssw for ATLAS, LHCb and CMS,
respectively). Even within a single framework, multiple different approaches
may exist which have been tailored for a specific use-cases.

Detailed detector performance studies typically need the most detailed
approach available, which comes at the cost of significant highly processing
times. Such studies can generally be performed to a sufficient level of
precision by using smaller samples sizes and/or a limited set of physics
processes. This is not the case for the production of Monte Carlo samples for
use in data analyses, where large samples covering a wide range of signal and
background processes are needed.

Using fully-detailed approaches for production of physics analysis Monte Carlo
samples can therefore sometimes become prohibitive, and this will increasingly
become the case for High-Luminosity LHC. Faster and more approximate
approaches therefore also need to be made available to be used in the
production of such samples, potentially with different approaches being used
for signal processes to those used for backgrounds, or the minimum bias
collision events used to represent the large number of additional ‘pile-up’
interactions per bunch crossing which accompany the signal process. In certain
circumstances, different levels of detail may be considered at different
positions within the detector, with a more accurate modeling applied for
crucial measurements close to the particle interaction point, or in areas of
high particle flux.

The faster, more approximate methods will typically rely on a parameterization
or templates derived from the fully-detailed approach, stand-alone
simulations, or from data. Appropriate and robust procedures for making
comparisons of key observables between the different approaches and also
crucially between experiment-specific and stand-alone software are therefore
highly beneficial for this process.

One aspect of silicon detector simulation which is becoming increasingly
important is the modeling of radiation damage. Different approaches to this
topic include111111See also this talk for an overview.:

  * •

A modified digitization approach has been employed by ATLAS, in which specific
adaptations are made within the code responsible for modeling charge
generation and transport within the sensor, in order to account for the
effects of radiation damage on the generated charges prior to their readout
RadDamageDigi_2019 .

  * •

CMS takes a different approach, in which _cluster templates_ are derived based
on comparing the properties of charge clusters produced in independent
simulations with and without radiation damage effects. These templates are
parameterized according to pertinent factors such as particle incident angle,
and position within the detector, in order that appropriate efficiency
corrections can be applied to channels within the clusters produced without
explicit radiation damage modeling applied.

  * •

LHCb uses an even faster approach whereby fundamental silicon parameters (such
as the diffusion) are tuned to data to reproduce the effective consequences of
radiation damage. For example, a multiplicative factor is used to reduce the
predicated collected charge and the diffusion length parameter is increased to
match data. Currently, these parameters are updated annually.

##  3 Challenges and Needs

###  3.1 Proprietary software and device properties

One of the main systematic uncertainties in device simulation is the doping
profile of a sensor. This can both be the doping profile of a sensitive layer
such as an epitaxial layer in a monolithic sensor, or the doping profile of an
implant munker_thesis . These properties of a sensor are often proprietary and
are adjusted in simulations through tuning the simulation to measurements.
Also, to understand the effect of the uncertainty, doping profiles of a sensor
varied in both electrostatic and transient simulations.

In addition, device simulation software like TCAD is proprietary, leaving the
possibility for use up to whether an institute can afford a license and no
possibility for writing extensions. This also creates a dependence on only one
correct installation of the software per institute, as personal installation
is not possible. The Weightfield2 weightfield2  package provides an open
source alternative to TCAD that requires much less computing power but is also
not a replacement for the quite accurate device simulation in TCAD.

###  3.2 Computational requirements of TCAD-simulations

Two-dimensional electrostatic and transient simulations with software like
TCAD can often provide a good approximation of the electric and weighting
fields in a three-dimensional sensor. A sensor design is translated into a
mesh of points, where computations are done for each point in the mesh. To get
a more accurate picture of a sensor, three-dimensional simulations are also
possible. These results, however, in a much larger amount of mesh points for
which each step of the finite element simulation is stored. To optimize the
number of mesh points, one can optimize the position of a representation mesh
of a sensor using boundary conditions where edge effects are avoided.

###  3.3 Combined TCAD and transient simulation tools

As described in 2.6, Monte Carlo simulation tools such as Allpix2 allpix2 ,
KDetSim kdetsim1 , and TCode Loi2020  can be combined with electrostatic and
transient simulations from device simulation tools like TCAD. Allpix2 has an
approximation of a linear electric field, but this is in the case for
monolithic sensors with small collection electrodes, for example, not enough
to describe the complex electric fields resulting from the sensor design. The
sensor design needs to be implemented precisely, and knowledge from planar
sensors cannot be transferred. A two-dimensional or three-dimensional electric
field as well as weighting field of a sensor can be fed into the Allpix2
software but here care needs to be taken for it to be in the correct format.
The Allpix2 package includes a converter for a certain mesh file format
produced with TCAD simulations and plans for an interface with Weightfield2
spannagel2021allpix2 , but an all-in-one package does not yet exist: there is
still a need to set up simulation chains.

###  3.4 Uncertainties in TCAD models

Uncertainties in TCAD models, both systematic and statistical, are not always
used and mentioned consistently, nor can they always be estimated well. As
described above, doping levels of silicon sensors can influence electric
fields, the depletion region, and signal response. These are, however, very
often proprietary and not known to the author of the simulation. Temperature
can influence leakage currents, and is also under certain circumstances not
known. Different charge carrier lifetimes can also influence leakage currents,
as well as signal response. The position of the grid on which the finite
element simulation is computed can, if not fine enough, yield inaccurate
results. Finally, using two or three dimensions affect the accuracy of the
simulation but also the computing time.

###  3.5 Unified radiation damage (TCAD) and annealing

The detailed mix of particle types and energies can have a strong impact on
the level of macroscopic radiation damage effects and even more on the damage
at microscopic level, in particular on the electric field profile inside the
sensor bulk Kramberger_2014 ; MikuzVertex2017 . So, even if a unified TCAD
radiation damage model would be desirable, still different models are needed,
depending on the type of irradiation. For example, from the experience of
silicon detectors at hadron colliders (LHC, HL-LHC) it would be good to derive
TCAD radiation damage models knowing in detail the mix of particle types and
energies that had crossed the detectors during the data taking periods. This
knowledge could be then used for the future experiments ad hadronic colliders
(FCC-hh), but once the expected cocktail of irradiation has been carefully
simulated.

At the moment no TCAD radiation damage model includes a parameterization for
annealing effects. Often the TCAD radiation damage models have been based on
data from irradiated sensor that underwent after irradiation some annealing
process before being measured Chiochia2005 ; FOLKESTAD201794 ; schwandt2019new
. The typical annealing process for irradiated test structures is the so-
called “optimal annealing” (80 minutes at 60∘ C), which assures achieving the
smallest depletion voltage for a fixed level of irradiation. In installed and
running detectors it is difficult to have a precise annealing scheme: if it is
true that in LHC experiments silicon trackers are operated cold ($<0^{\circ}$
C), the temperature during maintenance is more difficult to be controlled,
hence some unwanted annealing can take place. Simulating this effect in TCAD
turns out to be difficult as it could be already being taken into account by
the original model.

For what concerns experiments at hadronic colliders, in the simulation of
radiation damage to ATLAS Pixels RadDamageDigi_2019  an effective approach has
been adopted, which consisted in matching the predicted annealing effect moll-
thesis  on the effective doping concentration $N_{eff}$ with the concentration
of simulated deep-defects. The effect turned out not to be so large thanks to
the fact that ATLAS Pixels detector was kept cold most of the time during its
operation. This may not always be the case and for significant annealing
effects, more sophisticated approaches are required.

###  3.6 Improving the Hamburg Model Simulation

The Hamburg Model code is written in C++ and utilizes tools in ROOT. Both
leakage current and depletion voltage predictions using the Hamburg Model
require as inputs the temperature of the sensor, fluence, and duration of the
fluence application on the detector. The code iterates through all temporal
measurement points to predict the leakage current and depletion voltage at a
future time. The time taken for the iteration depends on the total duration of
the fluence application which is the operational period of the detector. This
can lead to unsupportable run-time for long-term predictions. However, the
code can be accelerated. The methods for accelerating the code are in
development.

Addition of charge collection efficiency prediction routines is also proposed.
Furthermore, geometric generalization that allows modeling of alternative
sensor geometries (such as 3D versus planar) is needed. A multi-threaded
implementation of the code will significantly improve the performance of the
code in iterating through all fluence points.

###  3.7 Damage Models

One central uncertainty that obfuscates the translation between predicted
fluence and sensor damage is the precision of sensor damage factors. These
factors describe how much damage is introduced by a particle of a given energy
and given type. Many testbeam measurements have been made to constrain these
parameters, but most measurements are not produced with uncertainties and many
factors at energy values away from beam energies are based on simulations or
other extrapolations. The RD50 Collaboration maintains a database of damage
factors121212https://rd50.web.cern.ch under NIEL., but a dedicated measurement
program is needed to modernize and improve them. See Ref. ATLAS:2021gld  for
additional challenges associated with sensor modeling and with current
unexplained systematic differences between data and predictions.

###  3.8 Challenges of TCAD radiation-defect modeling for future colliders

For more complete parameterizations of radiation-induced properties like
CCE($x$), test beam measured charge collection data of varying pitches of both
strip and pixel detectors as well as higher number of fluence points would be
required. Future efforts of the defect model developments should also include
further calibration with the indirect electric field information from the
measured edge-Transient Current Technique (TCT) data that enables the tuning
of the simulated $E$($z$) profile with increased precision.

The main tuning effort in the TCAD defect models discussed in Section 2.1 has
typically been focused to a fixed polarity and geometry device, leading to
compromised accuracy of the simulation results when these are varied. A more
general expression of device parameters would provide a starting point towards
unified defect models. Furthermore, effort should be put to tune the standard
parameter sets of the Sentaurus and Silvaco TCAD packages to produce
converging results of the basic sensor properties (that also include
uncertainties). This would lead to a more complete understanding and control
over the simulation process, and potentially to merging defect models between
the packages. It should not be taken for granted that basic silicon properties
are well-known ATL-INDET-PUB-2018-001 .

The essential objective of the combined bulk and surface damage simulations is
to stretch the validated defect modeling up to $\sim 2\times 10^{16}$
n${}_{\textrm{\small eq}}/$cm2 to account for the fluences of the pixel and
3D-columnar pixel detectors positioned closest to the vertex at the HL-LHC, as
well as for sensors utilized in other extreme-fluence future colliders (e.g.
$e^{+}$/$e^{-}$ and hadron colliders). Presently the validation limit of most
of the peer-reviewed published models is about $1\times 10^{15}$
n${}_{\textrm{\small eq}}/$cm2, with the exception of separate characteristics
like LGAD leakage current in Croci:2022ygd  and CCE of 300-$\upmu$m-thick
strip sensors and 3D-pixel sensors in Passeri2016  and Dawson:2764325 ,
respectively, in the $10^{16}$ n${}_{\textrm{\small eq}}/$cm2 fluence region.

All these TCAD-simulation R&D branches will require significant amount of
effort (time, budget, personnel, etc.) and will be difficult to realize
without a dedicated and well-organized collaboration.

## Acknowledgments

BN was supported by the Department of Energy, Office of Science under contract
number DE-AC02-05CH11231. XS was supported by the National Natural Science
Foundation of China (No. 1196114014).

##  Appendix A TCAD mathematical models

TCAD software suites are designed to simulate charge transport and electrical
polarization of semiconductor devices using a semi classical approach to
electronic transport. The essential set of equations, Fu2012 , takes into
account the Poisson equation, 7 for the electric field and charges densities,
the current continuity, eqs. 8, 9 and the trap occupation dynamics, eqs.10 to
12:

$\displaystyle\nabla\cdot(\epsilon\nabla\phi)$ $\displaystyle=$
$\displaystyle-q(p-n+N_{D}(1-f^{n}_{D})-N_{A}f^{n}_{A})-q\sum_{j}N_{tj}(\delta_{j}-f^{n}_{tj})$
(7) $\displaystyle\frac{\partial n}{\partial t}-N_{D}\frac{\partial
f^{n}_{D}}{\partial t}$ $\displaystyle=$
$\displaystyle\big{(}G_{net,n}-\sum_{j}R_{n}^{tj}-R_{au}\big{)}+\frac{1}{q}\nabla\cdot\overrightarrow{J_{n}}$
(8) $\displaystyle\frac{\partial p}{\partial t}+N_{A}\frac{\partial
f^{n}_{A}}{\partial t}$ $\displaystyle=$
$\displaystyle\big{(}G_{net,p}-\sum_{j}R_{p}^{tj}-R_{au}\big{)}-\frac{1}{q}\nabla\cdot\overrightarrow{J_{p}}$
(9) $\displaystyle N_{tj}\frac{\partial f^{n}_{tj}}{\partial t}$
$\displaystyle=$ $\displaystyle R_{n}^{tj}-R_{p}^{tj}\qquad\textrm{for each
trap j, where:}$ (10) $\displaystyle R_{n}^{tj}$ $\displaystyle=$
$\displaystyle c_{nj}nN_{tj}(1-f^{n}_{tj})-e_{nj}N_{tj}f^{n}_{tj}$ (11)
$\displaystyle R_{p}^{tj}$ $\displaystyle=$ $\displaystyle
c_{pj}pN_{tj}f^{n}_{tj}-e_{pj}N_{tj}(1-f^{n}_{tj})$ (12)

The summation term in Eq.(7) is $\rho_{trap}$, the trapped net charge, where
$\delta_{j}=1$ if the trap j is of donor type, $\delta_{j}=0$ if the trap j is
of acceptor type. In Eqs.(11, 12) $c_{nj},\,c_{pj}$ are the electron/hole
capture terms and $e_{nj},\,e_{pj}$ are the electron/hole emission terms for
trap $j$. $f^{n}_{tj}$ means the electron occupation fraction of trap $j$ and
$f^{n}_{D}$, $f^{n}_{A}$ are the occupancies for dopant donors and acceptors.
$G_{net,n},\,G_{net,p}$ are the net generation rates for electrons and holes,
including optical, radiation, impact ionization and other available generation
mechanisms. Correspondingly, $R_{au}$ is the Auger recombination and
$R_{n}^{tj},\,R_{p}^{tj}$ are the electron and hole recombination rates for
the trap $j$.

Different TCAD software providers offers slightly different set of equations
for the currents but all of them implement the drift-diffusion (DD) currents
for $\overrightarrow{J}_{n}$ and $\overrightarrow{J}_{p}$. DD is appropriate
for simulation of semiconductor particle detectors at constant temperature.
The semiconductor is not in equilibrium so it is necessary to define Quasi-
Fermi potentials, $\Phi_{n},\Phi_{p}$ grundmann . The DD currents definition
is, Chang1988 :

$\displaystyle\overrightarrow{J}_{n}$ $\displaystyle=$ $\displaystyle-
nq\mu_{n}\nabla\Phi_{n}=\mu_{n}(n\nabla E_{c}-\frac{3}{2}nKT\nabla
ln\,m_{n})+D_{n}(\nabla n-n\,\nabla ln\,\gamma_{n})$ (13)
$\displaystyle\overrightarrow{J}_{p}$ $\displaystyle=$ $\displaystyle-
pq\mu_{p}\nabla\Phi_{p}=\mu_{p}(n\nabla E_{v}-\frac{3}{2}nKT\nabla
ln\,m_{p})-D_{p}(\nabla p-p\,\nabla ln\,\gamma_{p})$ (14)

that considers the contributions from drift, diffusion and also spatial
variations of the effective mass and Fermi-Dirac statistics.

## References

  * (1) I. Dawson, _Radiation effects in the LHC experiments: Impact on detector performance and operation_.  CERN Yellow Reports: Monographs. CERN, Geneva, 2021, 10.23731/CYRM-2021-001.
  * (2) R. Radu et al., _Investigation of point and extended defects in electron irradiated silicon - Dependence on the particle energy_ , _J. Appl. Phys._ 117 (2015) 164503.
  * (3) V. Eremin, E. Verbitskaya and Z. Li, _The origin of double peak electric field distribution in heavily irradiated silicon detectors_ , _Nucl. Instr. $\&$ Meth. A_ 476 (2002) 556–564.
  * (4) V. Chiochia et al., _Simulation of heavily irradiated silicon pixel sensors and comparison with test beam measurements_ , _IEEE Trans. Nucl. Sci_ 52 (2005) 1067–1075.
  * (5) V. Eremin, E. Verbitskaya, A. Zabrodskii, Z. Li and J. Härkönen, _Avalanche effect in Si heavily irradiated detectors: Physical model and perspectives for application_ , _Nucl. Instr. $\&$ Meth. A_ 658 (2011) 145–151.
  * (6) E. Verbitskaya, V. Eremin and A. Zabrodskii, _Operational voltage of silicon heavily irradiated strip detectors utilizing avalanche multiplication effect_ , _JINST_ 7 (2012) C02061.
  * (7) R. Eber, _Investigations of new sensor designs and development of an effective radiation damage model for the simulation of highly irradiated silicon particle detectors_ , Ph.D. thesis, Karlsruhe Institute of Technology, 2013\.  http://digbib.ubka.uni-karlsruhe.de/volltexte/1000038403.
  * (8) M. Petasecca et al., _Numerical simulation of radiation damage effects in p-type silicon detectors_ , _Nucl. Instr. $\&$ Meth. A_ 563 (2006) 192–195.
  * (9) D. Pennicard et al., _Simulations of radiation-damaged 3D detectors for the Super-LHC_ , _Nucl. Instr. $\&$ Meth. A_ 592 (2008) 16–25.
  * (10) D. Passeri, F. Moscatelli, A. Morozzi and G. M. Bilei, _Modeling of radiation damage effects in silicon detectors at high fluences HL-LHC with Sentaurus TCAD_ , _Nucl. Instr. $\&$ Meth. A_ 824 (2016) 443–445.
  * (11) Å. Folkestad, K. Akiba, M. van Beuzekom, E. Buchanan, P. Collins, E. Dall’Occo et al., _Development of a silicon bulk radiation damage model for sentaurus tcad_ , _Nucl. Instr. $\&$ Meth. A_ 874 (2017) 94 – 102.
  * (12) J. Schwandt, E. Fretwurst, E. Garutti, R. Klanner, C. Scharf and G. Steinbrueck, _A new model for the TCAD simulation of the silicon damage by high fluence proton irradiation_ , 1904.10234.
  * (13) G. Verzellesi, G.-F. D. Betta and G. Pignatel, _Compact modeling of n-side interstrip resistance in p-stop and p-spray isolated double-sided silicon microstrip detectors_ , _IEEE Nucl. Sci. Symp._ 1 (2000) 25–27.
  * (14) C. Piemonte, _Device Simulations of Isolation Techniques for Silicon Microstrip Detectors Made on p-Type Substrates_ , _IEEE Trans. Nucl. Sci_ 53 (2006) 1694–1705.
  * (15) Y. Unno et al., _Optimization of surface structures in n-in-p silicon sensors using TCAD simulation_ , _Nucl. Instr. $\&$ Meth. A_ 636 (2011) 118–124.
  * (16) J. Zhang et al., _Study of radiation damage induced by 12 keV X-rays in MOS structures built on high-resistivity n-type silicon_ , _J. Synchrotron Rad._ 19 (2012) 340–346.
  * (17) T. Peltola, _Charge collection efficiency simulations of irradiated silicon strip detectors_ , _JINST_ 9 (2014) C12010.
  * (18) T. Peltola et al., _A method to simulate the observed surface properties of proton irradiated silicon strip sensors_ , _JINST_ 10 (2015) C04025.
  * (19) R. Dalal et al., _Combined effect of bulk and surface damage on strip insulation properties of proton irradiated n +-p silicon strip sensors_, _JINST_ 9 (2014) P04007.
  * (20) T. Peltola, _Simulation of radiation-induced defects_ , _PoS_ 031 (2015) (VERTEX2015).
  * (21) N. Akchurin et al., _Charge Collection and Electrical Characterization of Neutron Irradiated Silicon Pad Detectors for the CMS High Granularity Calorimeter_ , _JINST_ 15 (2020) P09031.
  * (22) F. Moscatelli, D. Passeri, A. Morozzi, S. Mattiazzo, G. F. Dalla Betta, M. Dragicevic et al., _Effects of Interface Donor Trap States on Isolation Properties of Detectors Operating at High-Luminosity LHC_ , _IEEE Trans. Nucl. Sci._ 64 (2017) 2259–2267.
  * (23) A. Morozzi, F. Moscatelli, D. Passeri and G. M. Bilei, _TCAD advanced radiation damage modeling in silicon detectors_ , _PoS_ 050 (2020) (VERTEX2019).
  * (24) A. Morozzi, F. Moscatelli, T. Croci and D. Passeri, _TCAD Modeling of Surface Radiation Damage Effects: A State-Of-The-Art Review_ , _Front. in Phys._ 9 (2021) 617322.
  * (25) T. Croci, A. Morozzi, F. Moscatelli, V. Sola, G. Borghi, G. Paternoster et al., _TCAD simulations of non-irradiated and irradiated low-gain avalanche diodes and comparison with measurements_ , _JINST_ 17 (2022) C01022.
  * (26) D. Massey, J. David and G. Rees, _Temperature dependence of impact ionization in submicrometer silicon devices_ , _Electron Devices, IEEE Transactions on_ 53 (10, 2006) 2328 – 2334.
  * (27) J. Calvo et al., _TRACS: a multi-thread transient current simulator for micro strips and pad detectors_ , _Nucl. Instr. $\&$ Meth. A_ 917 (2016) 77–85.
  * (28) F. Cenna et al., _Weightfield2: a fast simulator for silicon and diamond solid state detector_ , _Nucl. Instr. $\&$ Meth. A_ 796 (2015) 149–153.
  * (29) G.Kramberger, _Kdetsim - a ROOT based detector simulation package_ , in _27th RD50 Workshop_ , 2015, https://indico.cern.ch/event/456679/contributions/1126324/.
  * (30) W. Shockley, _Currents to conductors induced by a moving point charge_ , _J. Appl. Phys._ 9 (1938) 635–636.
  * (31) S. Ramo, _Currents induced by electron motion_ , _Proc. IRE_ 27(9) (2015) 584–585.
  * (32) T. Peltola et al., _Characterization of thin p-on-p radiation detectors with active edges_ , _Nucl. Instr. $\&$ Meth. A_ 813 (2016) 139–146.
  * (33) M. Alyari et al., _200 mm sensor development using bonded wafers_ , _JINST_ 16 (2021) T02002, [2006.04888].
  * (34) M. Moll, _Radiation damage in silicon particle detectors: Microscopic defects and macroscopic properties_ , Ph.D. thesis, Hamburg U., 1999.  http://www-library.desy.de/cgi-bin/showprep.pl?desy-thesis99-040.
  * (35) G. Lindstrom, _Radiation damage in silicon detectors_ , _Nucl. Instrum. Meth. A_ 512 (2003) 30–43.
  * (36) A. Chilingarov, _Temperature dependence of the current generated in si bulk_ , _JINST_ 8 (2013) P10003.
  * (37) S. Baranov, M. Bosman, I. Dawson, V. Hedberg, A. Nisati and M. Shupe, _Estimation of Radiation Background, Impact on Detectors, Activation and Shielding Optimization in ATLAS_ , tech. rep., CERN, Geneva, Jan, 2005.  http://cds.cern.ch/record/814823.
  * (38) G. Battistoni, S. Muraro, P. R. Sala, F. Cerutti, A. Ferrari, S. Roesler et al., _The FLUKA code: Description and benchmarking_ , _AIP Conf. Proc._ 896 (2007) 31–49.
  * (39) GEANT4 collaboration, S. Agostinelli et al., _GEANT4–a simulation toolkit_ , _Nucl. Instrum. Meth. A_ 506 (2003) 250–303.
  * (40) J. Allison et al., _Geant4 developments and applications_ , _IEEE Trans. Nucl. Sci._ 53 (2006) 270.
  * (41) J. Allison et al., _Recent developments in Geant4_ , _Nucl. Instrum. Meth. A_ 835 (2016) 186–225.
  * (42) ATLAS collaboration, G. Aad et al., _Measurements of sensor radiation damage in the ATLAS inner detector using leakage currents_ , _JINST_ 16 (2021) P08025, [2106.09287].
  * (43) M. Alnaes et al., _The FEniCS Project Version 1.5_ , _Archive of Numerical Software_ 3 (2015) .
  * (44) F. Cenna, N. Cartiglia, M. Friedl, B. Kolbinger, H.-W. Sadrozinski, A. Seiden et al., _Weightfield2: A fast simulator for silicon and diamond solid state detector_ , _Nucl. Instr. $\&$ Meth. A_ 796 (2015) 149–153.
  * (45) T. Yang et al., _Time Resolution of the 4H-SiC PIN Detector_ , _Frontiers in Physics_ 10 (2022) .
  * (46) Y. Tan et al., _Timing Performance Simulation for 3D 4H-SiC Detector_ , _Micromachines_ 13 (2022) .
  * (47) R. J. Lipton, _LGAD Single Event Burnout Studies_ , tech. rep., FNAL, United States, Aug, 2021.  https://doi.org/10.2172/1841397.
  * (48) D. Dannheim, K. Dort, D. Hynds, M. Munker, A. Nürnberg, W. Snoeys et al., _Combining TCAD and Monte Carlo methods to simulate CMOS pixel sensors with a small collection electrode using the Allpix 2 framework_, _Nucl. Instr. $\&$ Meth. A_ 964 (2020) 163784.
  * (49) M. Swartz, V. Chiochia, Y. Allkofer, D. Bortoletto, L. Cremaldi, S. Cucciarelli et al., _Observation, modeling, and temperature dependence of doubly peaked electric fields in irradiated silicon pixel sensors_ , _Nucl. Instr. $\&$ Meth. A_ 565 (2006) 212–220.
  * (50) M. Swartz, _A Detailed Simulation of the CMS Pixel Sensor_ , tech. rep., CERN, Geneva, Jul, 2002.  https://cds.cern.ch/record/687440.
  * (51) M. Swartz, D. Fehling, G. Giurgiu, P. Maksimovic and V. Chiochia, _A new technique for the reconstruction, validation, and simulation of hits in the CMS pixel detector_ , _PoS_ 035 (2007) (VERTEX2007).
  * (52) V. Chiochia, M. Swartz, D. Fehling, G. Giurgiu and P. Maksimovic, _A novel technique for the reconstruction and simulation of hits in the cms pixel detector_ , in _2008 IEEE Nuclear Science Symposium Conference Record_ , pp. 1909–1912, 2008, DOI.
  * (53) CMS collaboration, S. Chatrchyan et al., _Description and performance of track and primary-vertex reconstruction with the CMS tracker_ , _JINST_ 9 (2014) P10009, [1405.6569].
  * (54) M. Benoit, J. Idarraga and S. Arfaoui, _The AllPix Simulation Framework: Generic Geant4 implementation for pixel detectors_ , _Last accessed 2022-01-25._ https://twiki.cern.ch/twiki/bin/view/Main/AllPix#Latest_Version (2017) .
  * (55) S. Spannagel, K. Wolters, D. Hynds, N. Alipour Tehrani, M. Benoit, D. Dannheim et al., _Allpix2: A modular simulation framework for silicon detectors_ , _Nucl. Instr. $\&$ Meth. A_ 901 (2018) 164–172.
  * (56) S. Spannagel and P. Schütze, _Allpix 2 – silicon detector monte carlo simulations for particle physics and beyond_, 2112.08642.
  * (57) M. Aaboud et al., _Modelling radiation damage to pixel sensors in the ATLAS detector_ , _JINST_ 14 (2019) P06012.
  * (58) M. Munker, _Test beam and simulation studies on High Resistivity CMOS pixel sensors_ , Ph.D. thesis, Bonn U., May, 2018.  https://cds.cern.ch/record/2644054.
  * (59) A. Loi and A. Contu, _Tcode: A new multithread simulator for silicon sensors in hep applications_ , in _2019 IEEE Nucl. Sci. Symp. (NSS/MIC)_ , pp. 1–4, 2019, DOI.
  * (60) G. Kramberger, V. Cindro, I. Mandić, M. Mikuž, M. Milovanović and M. Zavrtanik, _Modeling of electric field in silicon micro-strip detectors irradiated with neutrons and pions_ , _JINST_ 9 (2014) P10016.
  * (61) M. Mikuž, G. Kramberger, V. Cindro, I. Mandić and M. Zavrtanik, _Extreme radiation tolerant sensor technologies_ , 2017.  https://indico.ijs.si/event/875/contributions/2657/attachments/946/1201/Extreme-Vertex-Sep17.pdf.
  * (62) ATLAS Collaboration, _Modeling the Mobility and Lorentz angle for the ATLAS Pixel Detector_ , _ATL-INDET-PUB-2018-001_ (2018) .
  * (63) S. Li and Y. Fu, _3D TCAD Simulation for Semiconductor Processes, Devices and Optoelectronics_.  Springer Science+Business Media. Springer, New York,, 2012, https://doi.org/10.1007/978-1-4614-0481-1.
  * (64) M. Grundmann, _The Physics of Semiconductors_.  Graduate Texts in Physics. Springer, New York, 2006, https://doi.org/10.1007/978-3-030-51569-0.
  * (65) K.M.Chang, _A consistent model for carrier transport in heavily doped semiconductor devices_ , _Semiconductor Science and Technology_ 3(8) (1988) 766–772.