{'fraction_non_alphanumeric': 0.02700416152192802, 'fraction_numerical': 0.011455618054602826, 'mean_word_length': 4.964948172394981, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulations of Silicon Radiation Detectors for High Energy Physics Experiments

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.


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. No one model 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.

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. [1].
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 degrade various operational properties of a silicon sensor both inside the Si-bulk and at the surface on the Si/SiO 2 -interface 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 f ) with/without 1-2 acceptor and donor type interface-trap (N 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-detectorsystem 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 [1], 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 (∼ 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 radia-tion 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 × 10 16 n eq /cm 2 , expected for future colliders, are considered in Sections 3.5, 3.6 and 3.8.

TCAD Simulations for Detector Properties
Radiation levels above ∼ 10 13 n eq /cm 2 introduce observable damage to the crystal structure of a silicon sensor. Fluences beyond 1 × 10 14 n eq /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 SiO 2 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 [2], 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 software 1 frameworks Synopsys Sentaurus™ 2 and Silvaco Atlas™ 3 .
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 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 [3][4][5][6][7], three-level [8][9][10][11], as well as a fivelevel models [12]. Simulating the effects of the surface damage on Si/SiO 2 -interface irradiated by charged hadrons, gammas, or X-rays -like the change of inter-electrode resistance (R int ) and capacitance (C 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 f ) [13][14][15] or by three interface traps (N it ) with parameters matching the measurements of X-ray irradiated MOS-capacitors [16]. Results from a further model that combines the two approaches is presented in Figure 2, where the characteristics of the CV -curves of a γ-irradiated MOS-capacitor are only reproduced by simulation by including in addition to N f both acceptor and donor-type deep N it (N it,acc,don ) at the Si/SiO 2 -interface ( Figure 2a). By using both N f and N it,acc,don as an input for an R int -simulation with typical CMS High Granularity Calorimeter (HG-CAL) isolation implant parameters in Figure 2b, the experimentally observed R int -values are quantitatively reproduced (atoll and common p-stops with STD N ps ). By excluding N it,acc,don from the simulation, the pads become either shorted (STD N ps , N 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 × STD N ps , N 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 it with 2 µm depth distribution from the surface (Sentaurus [17,18]), and another approach where two bulk defect levels are complemented by two acceptor-type N it (Silvaco [19]). More details on the development of the defect models are described in e.g., [20].  Figure 1. Examples of bulk-damage modeling [21]. (a) Measured and simulated IR-laser induced transient signals for 300-µm-thick p-on-n pad-sensor (300N) after neutron irradiation to the fluence of (6.1 ± 0.5) × 10 14 n eq /cm 2 . The simulation applied neutron defect model [7] with a fluence of 6.0 × 10 14 n eq /cm 2 . The sensor parameters used in the simulation were extracted from CV /IVmeasurements. (b) Measured and simulated evolution of charge collection efficiency (CCE) with voltage and 1-MeV equivalent neutron fluence for 300-µm thick sensors at -30 • C. Simulated results use dotted, dashed and solid curves for 600 V, 800 V and 1 kV, respectively.
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 SiO 2 , 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 [22].
In the original Perugia model [10,22], 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 refereed to as the Perugia surface model 2019 [23,24]. The surface damage effects can be mainly described by three parameters: the oxide charge density (N ox ) and the two interface trap states (N it ) for donors and acceptors. The values of the above quantities can be extracted from highfrequency 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 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) (  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) [25]. 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 [26].
This last version of the Perugia model is able to reproduce the behavior of 3D-devices as well.

TCAD Alternatives
There are simpler software alternatives to TCAD like TRACS, [27], Weightfield2, [28] or KDetSim 4 [29]. All of them solve the Poisson equation (A.1) for a charge distribution N eff . Instead of calculating N 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, [30,31]: 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.

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. 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 [33].

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 the Hamburg Model [34]. A linear relationship between particle fluence and leakage current applies to silicon sensors. The effective fluence, Φ eq , is the number of 1-MeV neutrons applied to a sensor of surface area 1 cm 2 that cause damage equivalent to that of all particles that went through the sensor. The linear relationship between effective fluence and ∆I, the change in leakage current at fluence Φ eq relative to the value before irradiation, is given as where V is the depleted volume of the silicon sensor and α is the current-related damage coefficient [35] 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 eff [36]. 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 [37,38] and Geant4 [39][40][41], as given by the following equation: where E eff = 1.21 eV is used for the effective silicon band gap energy, k B is the Boltzmann constant, and T R = 21 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 [42] 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.

Modeling Depletion Voltage and Charge Collection Efficiency
The full depletion voltage V dep of a sensor depends linearly on the doping concentration N 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 where q is the charge of the electron, is the dielectric constant, and 0 is the vacuum permittivity. The doping concentration is affected by the fluence, and the change in doping concentration ∆N eff can be expressed as where N A , N Y , and N C are the changes in the doping concentration due to beneficial annealing, reverse annealing, and stable damage respectively. The Hamburg Model 5 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. 2.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 [34] and are not suitable for predictions of depletion voltage of detectors designed for future conditions. As an illustration of this limitation, a comparison [1] of depletion voltage prediction to measured depletion voltage of the ATLAS Insertable B-Layer is shown in Figure 9. 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.

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) 6 has been developed as a python module. In RASER, the electric field and weighting potential are calculated by FEniCS [43] -an open-source computing platform for solving partial differential equations. The simulation of electrons from 90 Sr 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 [30,31]. The readout electronics used a simplified current amplifier as in [44]. 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 × 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 ± 1) ps and (73 ± 1) ps, respectively [45].
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 × 1 cm, and the thickness is 350 µm, and the effective concentration of n-type substrate is set to 1 × 10 13 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 90 Sr 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 keep constant even with temperature 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 [46]. Future development of RASER will include the LGAD-SiC as well as the pixel-timing performance for SiC devices.

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 systematic of these effects, including linear energy transfer threshold and the effect of the field strength on breakdown margin [47]. Figure 14 shows transient curves of temperature and cathode current for a 50-µm-thick diode operated at 625 V struck by particles with ionization between 32 and 140 MIPs.  Figure 14. Transient simulation of junction temperature (a) and cathode current (b) for a 50µ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.

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 [48].
One example of a software package for drift and diffusion modeling of carriers in silicon sensors is PixelAV [4,49,50], 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 [51][52][53]. In the templates created with PixelAV provide both corrections to the hit position and the cluster charge.
Packages similar to PixelAV are Allpix [54] 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 Allpix 2 [55]. 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. The Allpix 2 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. The packages now also include different charge carrier mobility models and the possibility to load doping profiles for doping-dependent charge carrier lifetime calculations [56], and foresees to include trapping of charges.

Full detector systems
Within large HEP collaborations, such as those of the LHC-experiments, silicon detector simulation is typically performed using dedicated, proprietary 7 software written within the experiment's software framework (for example Athena 8 , Gauss 9 , and CMSSW 10 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 include 11 : • 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 [57].
• 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.

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 [58]. 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 [44] 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.

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. This 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.

Combined TCAD and transient simulation tools
As described in 2.6, Monte Carlo simulation tools such as Allpix 2 [55], KDetSim [29], and TCode [59] can be combined with electrostatic and transient simulations from device simulation tools like TCAD. Allpix 2 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 Allpix 2 software but here care needs to be taken for it to be in the correct format. The Allpix 2 package includes a converter for a certain mesh file format produced with TCAD simulations and plans for an interface with Weightfield2 [56], but an all-in-one package does not yet exist: there is still a need to set up simulation chains.

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.

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 [60,61]. 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 [4,11,12]. The typical annealing process for irradiated test structures is the socalled "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 • 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 [57] an effective approach has been adopted, which consisted in matching the predicted annealing effect [34] on the effective doping concentration N ef f 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.

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.

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 factors 12 , but a dedicated measurement program is needed to modernize and improve them. See Ref. [42] for additional challenges associated with sensor modeling and with current unexplained systematic differences between data and predictions.

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 [62].
The essential objective of the combined bulk and surface damage simulations is to stretch the validated defect modeling up to ∼ 2 × 10 16 n eq /cm 2 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×10 15 n eq /cm 2 , with the exception of separate characteristics like LGAD leakage current in [25] and CCE of 300-µm-thick strip sensors and 3D-pixel sensors in [10] and [1], respectively, in the 10 16 n eq /cm 2 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. 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, [63], takes into account the Poisson equation, A.1 for the electric field and charges densities, the current continuity, eqs. A.2, A.3 and the trap occupation dynamics, eqs.A.4 to A.6: 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 tj n , R tj p 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 implements the drift-diffusion (DD) currents for − → J n and − → 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, Φ n , Φ p [64]. The DD currents definition is, [65]: nKT ∇ln m n ) + D n (∇n − n ∇ln γ n ) (A.7) that considers the contributions from drift, diffusion and also spatial variations of the effective mass and Fermi-Dirac statistics.