Development of a silicon bulk radiation damage model for Sentaurus TCAD

This article presents a new bulk radiation damage model for p-type silicon for use in Synopsys Sentaurus TCAD. The model is shown to provide agreement between experiment and simulation for the voltage dependence of the leakage current and the charge collection efficiency, for fluences up to 8 × 1015 1 MeV neq∕cm. © 2017 CERN for the benefit of the Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
In particle physics experiments, silicon detectors are often operated in harsh radiation environments, and understanding the impact of radiation damage on the detector performance is key to their successful operation. Device simulations using Technology Computer Aided Design (TCAD) software packages are useful tools for investigating the effects of radiation damage, in particular for linking macroscopic observables to what is happening on small scales inside the sensor bulk.
In the following a radiation damage model for -type silicon, implemented in Synopsys Sentaurus Device 2 TCAD, is presented. The model has been developed in the context of the R&D programme for the upgraded LHCb Vertex Locator (VELO), which will be installed in the LHCb experiment at CERN in 2019/2020 [1]. The model aims to reproduce charge collection efficiencies (CCE) and current-voltage (I-V) curves up to a fluence of 8 × 10 15 1 MeV n eq /cm 2 , the expected maximum fluence after an integrated luminosity of 50 fb −1 . The model is validated using measurements on irradiated -on-pixel sensors from Hamamatsu Photonics K.K. 3 These sensors have a thickness of 200 μm and a pixel cell size of 55 × 55 μm 2 , and feature -stop isolation between pixels.
Radiation damage models for Synopsys Sentaurus TCAD, of varying scope, have been developed in the past by different groups [2][3][4][5][6]. Differences between the present model and other models with similar range of validity in terms of fluence, in particular the Perugia model [3], are discussed in Section 2.4.

Simulations
The Sentaurus Device program allows for solving the Poisson and carrier continuity equations on two-dimensional (2D) and threedimensional (3D) discretised semiconductor structures using finite element methods. In this work, two types of simulations were performed: • steady-state simulations, where leakage current and electric field as function of voltage were simulated by solving the stationary Poisson and charge transport equations, • transient simulations, where the time dependent Poisson and charge transport equations are solved for a given initial charge distribution. Fig. 1. Close-up of the pixel region of a (left) 2D geometry with three pixels and (right) a 3D geometry with four quarter pixels. The 2D mesh used in CCE simulations contains two additional pixels, i.e. a total of five. The sensors simulated in this work consist of a boron-doped bulk in ⟨100⟩ orientation with phosphorus-doped 39 μm wide pixel implants and -stop regions (with a higher concentration of boron) between the implants. A layer of SiO 2 is placed on top of the bulk, and a positive surface charge with a density of 1 × 10 12 cm −2 is applied on the Si-SiO 2 interfaces. Measurements reported in the literature show that the oxide charge saturates at around this value [7]. The doping profile is not known in detail from the manufacturer. The bulk doping concentration used in the simulation was tuned to reproduce the depletion voltage of the sensors to which the simulations are compared, = 140 V. For the other parameters of the doping profile, which are less critical for simulating bulk leakage current and CCE, order of magnitude estimates were used. The doping profile is visualised in Fig. 1 and the main parameters are listed in Table 1.

Meshing
The number of grid points was chosen as a compromise between computing time and precision. Finer meshes were tried in several cases without significant changes in the results. Table 1 lists the number of mesh points.

Symmetries and boundary conditions
Voltage boundary conditions are imposed on the electrode-silicon and electrode-oxide interfaces, whereas Neumann boundary conditions are applied at all other mesh boundaries.
A stationary solution of the Poisson and drift-diffusion equations will reflect the periodicity of the pixel matrix. For simulating the bulk leakage current, it is therefore sufficient to simulate a single pixel and scale the result to the active area of the sensor. In order to simulate the leakage current of a sensor after exposure to the non-uniform irradiation profile discussed in Section 3.2, separate simulations were performed for every fluence level; the resulting leakage currents were then scaled to the surface area that had been exposed to the respective fluence, and the I-V curves of all regions were summed.
For 2D simulations of the collected charge, a mesh containing five pixels was used, with the ionising particle passing through the central pixel. Fig. 1 (left) shows three of the pixels in the five pixel mesh. For highly irradiated detectors, it is important to include neighbouring pixels since charge trapping causes a non-negligible net charge to be induced on neighbouring pixels. Five pixels were found to be sufficient to account for this effect.
The geometry used for charge collection simulations in 3D is illustrated in Fig. 2, which shows a top view of a 3 ×3 pixel region with an ionising particle going through the central pixel. The dark square intersecting a quarter of the deposited charge shows the area implemented in the 3D mesh. The particle goes through the corner of the mesh, and the deposited charge is one fourth of the 2D case. For reasons of symmetry, the solutions in each quadrant should be identical except for numerical effects. Instead of simulating the larger area with a dashed boundary, it is therefore sufficient to simulate one quadrant (with Neumann boundary conditions) and scale all currents by a factor of four. It should be noted that this argument only applies to tracks passing through the pixel centre, which is the case for all simulations discussed below. In the 3D mesh considered, only quarter neighbouring pixels are included-in contrast to two full neighbouring pixels on each side in 2D. The difference between 2D and 3D simulation results was used for assigning a systematic uncertainty to the 2D CCE simulations.

Primary ionisation
The ionisation pattern produced by a charged particle crossing the detector is modelled in terms of a cylindrically symmetric, continuous  Table 3 Sensors used in this work. For uniformly irradiated sensors the value in the third column is simply the fluence, while for non-uniform profiles it corresponds to the fluence in the area 0 mm < < 5 mm. charge distribution with a constant density (80 electron-hole pairs per μm) along the direction of the track, and a Gaussian distribution (1 μm standard deviation) in the transverse direction. The collected charge is given by the integrated current (after subtracting the leakage current) on all pixels that cross a threshold of 1000 electrons (the threshold used in data taking with the tested sensors). The integration time is 25 ns; integrating for a longer time was found to make only a negligible difference.

Sensor
Only perpendicular tracks passing through the pixel centre were simulated in both 2D and 3D. For tracks passing through the interpixel region, the charge collection properties are more sensitive to the modelling of surface damage, which is outside the scope of this work.

Physics models
In this work, the drift-diffusion model has been used, which implies that the temperature of the whole device remains constant. The continuity equations contain one source term for every defect level, and an additional source term for avalanche generation. The defect source terms are given by the Shockley-Read-Hall generation-recombination expression [8,9]. Sentaurus TCAD allows for the use of neutral trap levels for current generation, but these have not been used. Other physics models taken into account include Fermi-statistics, avalanche multiplication (Van Overstraeten-De Man model [10]), band gap narrowing (Slotboom model [11]), high field mobility saturation and trap assisted tunnelling (Hurkx model [12]). Detailed descriptions of these models can be found in the Sentaurus Device User Guide [13] and the references therein.

Radiation damage modelling
Developing a TCAD radiation damage model consists in defining a set of defect states, characterised by their location (energy level) in the band gap, electron and hole capture cross-sections ( , ℎ ), concentration and type (i.e. whether they are a donor or an acceptor). In theory, one could implement all defect levels that have been measured experimentally, but this approach is at present computationally prohibitive. Alternatively, one can define a few effective defect states and tune their parameters so that the model reproduces experimental observations. In this work the latter approach is used.
The two energy levels (defects 1 and 2 in Table 2) proposed by Eremin et al. [14] were used as a starting point. These levels, sometimes called the EVL levels, comprise one donor and one acceptor and are known to reproduce the double junction electric field effect [14][15][16]. Eber has further shown that agreement with measured I-V curves, and to some extent CCE (up to 1 × 10 15 1 MeV n eq /cm 2 ), can be achieved by using only these two energy levels [5]. These levels has also been combined with surface defects to model surface effects by Peltola et al. [6].
In addition to the EVL levels, a third defect was introduced. The procedure used to tune these three defects is outlined below.
• The defect state concentrations are assumed to scale linearly with 1 MeV neutron equivalent fluence with a proportionality factor (introduction rate) . • One of the irradiated sensors (assembly S4 in Table 3, uniformly irradiated to a fluence of 4 × 10 15 1 MeV n eq /cm 2 ), was selected as a reference. • The cross-sections and introduction rates of the two EVL levels were tuned to reproduce the measured I-V curve of the reference sensor. • A second acceptor close to the conduction band, corresponding roughly to the position of the A-centre defect state [17], was then introduced. This ''shallow'' acceptor (only 0.2 eV from ) has only a minor influence on the current generation and space charge, so that it allows for tuning the CCE independently of the behaviour of the I-V curves.
• The parameters of the second acceptor (defect 3 in Table 2) were tuned so that at one given voltage the simulated CCE agrees with the measured CCE of the reference sensor. Varying the hole capture cross-section ℎ within reasonable limits has a negligible effect on the CCE since the probability of hole capture is already very low due to the large distance from the valence band. To limit the number of degrees of freedom when scanning the parameter space, ℎ was chosen to have the same value as the electron capture cross-section . In addition, as a check, a least square fit of the simulated CCE curve with respect to the measured CCE curve of the reference sensor was performed using the introduction rate of the shallow acceptor as a fit parameter. In doing this, the introduction rate came out only 2% higher than with the method of tuning at one voltage. Table 2 summarises the parameters of the model used in this work. The cross-sections of the deep defects can be seen to be larger than the values used by Eremin et al. [14] (1 × 10 −15 cm −2 ).
The most recent Perugia model [3] aims to be valid up to fluences of 2 × 10 16 1 MeV n eq /cm 2 and is a natural basis for comparison. Both the model presented here and the Perugia model contain three bulk defect levels and are tuned for -type silicon. The model used in this paper differs from the Perugia model in that it aims only to reproduce bulk effects, while the Perugia model also includes surface effects. Furthermore, our proposed model is compared to different types of measurements, namely the voltage dependence of both the current and charge collection efficiency. It is furthermore based on the trap levels from the EVL-model that are also used by Eber and CMS, which differs from the traps used in the Perugia model. Both models use two acceptors and one donor, but their parameters are different. While the Perugia model contains two deep acceptors, our model contains one shallow and one deep acceptor.

Sensitivity analysis
The parameters of the ''deep'' defects (i.e. the traps with energy levels close to the middle of the band gap) are highly correlated and the effects of different trap states are not simply additive. In order to estimate the sensitivity of the CCE and I-V curves to uncertainties in the   and approximately causes an offset in the CCE curve. Changing ℎ has no effect and was set equal to in the tuning procedure.

Irradiations
Both sensors irradiated with neutrons and sensors irradiated with protons have been used in this work. Proton-irradiated sensors were irradiated at the KIT cyclotron 4 with ∼23 MeV protons. Most of the proton-irradiated sensors were exposed to a non-uniform, piece-wise constant fluence profile intended to mimic the profile expected in the upgraded VELO. The fluence increases in seven steps up to either 4 × 10 15 1 MeV n eq /cm 2 or 8 × 10 15 1 MeV n eq /cm 2 , depending on the sample. The uncertainty in the fluence calibration is estimated to be ±20% [19]. Fig. 6 shows the profile for a sensor with a maximum fluence of 4 × 10 15 1 MeV n eq /cm 2 . For the sensor with a maximum fluence of 8 × 10 15 1 MeV n eq /cm 2 , all values are doubled.
The neutron irradiated sensor were irradiated at JSI Ljubljana 5 to a uniform fluence of 8 × 10 15 1 MeV n eq /cm 2 . The uncertainty in the fluence has been estimated to be ±10% [20].
The annealing times of the proton sensors are estimated to between one and two days in room temperature, while the neutron sensor is estimated to have annealed about five days in room temperature.
A list of tested sensor assemblies and their irradiation profiles is shown in Table 3.

Leakage current measurements
The laboratory setup used for measuring the I-V curves comprises a vacuum tank where the sensors were mounted on a Peltier-cooled copper support plate. The heat generated by the Peltier cooler was removed using a chiller circulating a water/glycol fluid underneath the hot side of the Peltier. The temperature was monitored by a PT100 sensor glued on the copper block. In order to determine a relation between the temperature of the copper block and the temperature of the silicon itself, a PT100 was glued on top of a test sensor. A linear relationship between sensor temperature and copper temperature was found, with the temperature difference increasing towards lower temperatures.

Charge collection efficiency measurements
The CCE was measured in test beams at the H8 secondary beam line of the CERN Super Proton Synchrotron (SPS), using a beam of charged hadrons with a momentum of 180 GeV / . The device under test was placed in the centre of the Timepix3 telescope and oriented perpendicularly with respect to the beam. In order to be able to compare the measured charge collection efficiency with the simulated scenario discussed in Section 2.2, only clusters associated to telescope tracks passing through a window of ±5 μm around the pixel centre were used in the analysis. In addition, the arrival time of the clusters was required to be within ±200 ns of the time of the corresponding telescope track. 5 http://www.rcp.ijs.si/ric/description-a.html. To convert the time-over-threshold (ToT) measurements provided by the Timepix3 ASIC to collected charge, calibration curves obtained from testpulses was used. These curves were determined for each pixel by injecting delta pulses of known amplitude 100 times into the ASIC frontend, and fitting the average ToT as a function of the injected charge to a function where 0 , 1 , , are fit parameters. A typical example of a ToT calibration curve for a single pixel is shown in Fig. 7. The uncertainty in the testpulse calibration constitutes the dominating contribution to the systematic error of the charge collection measurements. The error was determined by comparing the testpulse calibration curves to calibration measurements using a 241 Am source [21], and was found to be ≲5%. The most probable value (MPV) of the collected charge was found by fitting the charge spectrum with a Landau distribution convoluted with a Gaussian. One of the sensors discussed in this work (assembly S6 in Table 3) was characterised in beam tests before and after irradiation, and the pre-irradiation MPV at overdepletion, 15.8× 10 3 − , was used for normalising the CCE of all irradiated sensors.

Leakage current
The simulated I-V curves presented below were obtained using a 2D geometry, as the results from 2D and 3D simulations were found to be virtually identical. The measured I-V curve used for tuning the model ( = 4 × 10 15 1 MeV n eq /cm 2 , = −31.1 • C) is shown in Fig. 8, together with the simulated I-V curve. Simulation and experiment agree (except at low bias voltage), which is expected since this curve was used for tuning the model. Fig. 9 compares predicted I-V curves at = 8 × 10 15 1 MeV n eq /cm 2 with measurements from a neutron-irradiated sensor, at two different temperatures. The simulated curves shows more structure above ∼200 V than the measured ones, but the magnitudes of the current and the approximate slopes are in agreement. The level of agreement is better at −31.8 • C (which is approximately where the model was tuned) than at −37.9 • C, where the simulated current at 1000 V bias is ∼25% lower than the measured one. Fig. 10 compares measurements and predictions for non-uniformly irradiated sensors, at temperatures ∼6 − 7 • C lower than the temperature at which the model was tuned. The measurements for the three sensors with a maximum fluence of 4 × 10 15 1 MeV n eq /cm 2 agree with simulation. For the sensor irradiated to a maximum fluence of 8 × 10 15 1 MeV n eq /cm 2 , the simulated curve exhibits some features that are not present in the measurement, but the order of magnitude of the current is reproduced. The breakdown visible in the measured I-V curves (a) = −31.8 • C.
(b) = −37.9 • C.    is not expected to be reproduced since surface damage and edge effects were not included in the simulation. As can be seen from Fig. 11, the simulated I-V curve shows a transition from saturating to non-saturating behaviour as the fluence increases. For the I-V curves that do exhibit saturating behaviour, the damage parameter can be calculated, where = ( ) − ( = 0) is the increase in leakage current after irradiation and is the depleted volume (which is identical to the physical bulk volume in the low-fluence regime). Fig. 12 shows the simulated saturation current per volume as a function of fluence at −1000 V. The simulation was performed at = −38.1 • C, but for comparison with literature data the simulated leakage current was scaled to = 21 • C, using the expression for the temperature dependence of the bulk current [22], where ,eff = 1.214 eV [23] is the effective band gap after irradiation and is the Boltzmann constant. The damage parameter at 21 • C is found to be sim = (4.329 ± 0.001) × 10 −17 A cm −1 . Literature values for at 21 • C after 0.3 and ten days of annealing at room temperature are (6.40±0.43) × 10 −17 A cm −1 and (4.32±0.29) × 10 −17 A cm −1 , respectively [24]. The damage parameter predicted by the model can thus be considered in agreement with what is expected for a sensor after around ten days of annealing at room-temperature.  Fig. 13 shows that for estimating the collected charge at = 4 × 10 15 1 MeV n eq /cm 2 a 2D strip-like model represents a close approximation of the full 3D pixel geometry, given that one is considering a particle passing through the centre of a pixel. The deviation is largest at 400 V where ∼380 more electrons are collected in 2D, while the average deviation is ∼260 − . The 2D model can be seen to consistently collect more charge, which may be attributed to the fact that fewer neighbouring pixels are included in the 3D model. It is generally seen that including nearest and next to nearest neighbours in 2D will cause an increase in collected charge on the central pixel and thus the CCE. In 3D, the neighbouring pixels are only partly included, and next to nearest neighbours are not included at all.

Collected charge
The measured charge collection efficiency for the reference sensor ( = 4 × 10 15 1 MeV n eq /cm 2 ) used for tuning the model is shown in Fig. 14(a) together with the simulated CCE (at = −32.1 • C). As discussed in Section 2.4, only one CCE measurement (the point at 700 V bias) was used for tuning the model; the voltage dependence is predicted by the model and is found to be in agreement with the data. Fig. 14(b) shows simulated and measured CCE as a function of bias voltage at a fluence of 8 × 10 15 1 MeV n eq /cm 2 . The measured CCE is close to the upper end of the fluence error band of the simulation (i.e. the simulated CCE for a fluence at the lower bound of the estimated uncertainty). It should be noted that in Fig. 14(a), the simulated CCE is compared to a proton-irradiated sensor, whereas at 8 × 10 15 1 MeV n eq /cm 2 it is compared to a neutron-irradiated sensor.
The temperature of the sensors in test beam was known to within ±3 • C, and the variation in the simulated CCE resulting from varying the temperature within this range was factored into the simulation error bars. At 4 × 10 15 1 MeV n eq /cm 2 , the contribution of the temperature uncertainty to the error bars is negligible, while at 8 × 10 15 1 MeV n eq /cm 2 it becomes significant at bias voltages > 500 V. In other words, the temperature dependence of the simulated CCE increases with increasing fluence and bias voltage. An explanation for this observation is given in Section 5.

Discussion
Fig . 15 shows the absolute value of the electric field in the pixel centre as a function of the depth in the sensor for fluences between 1 × 10 15 1 MeV n eq /cm 2 and 8 × 10 15 1 MeV n eq /cm 2 . Comparing Figs. 11 and 15, one can see that the double junction develops in the same fluence range where the leakage current exhibits a transition from saturating to non-saturating behaviour. As can be seen from Fig. 15, at the fluence levels characterised by a double junction field, there are regions in the sensor where the electric field exceeds 1 × 10 5 V cm −1 , which is the order of magnitude at which avalanche multiplication starts (a) 4 × 10 15 1 MeV n eq /cm 2 .

Fig. 14.
Comparison of measured and simulated CCE as a function voltage at = 4 × 10 15 1 MeV n eq /cm 2 (proton irradiation) and = 8 × 10 15 1 MeV n eq /cm 2 (neutron irradiation) using a 2D mesh. Only the data point at 700 V and 4 × 10 15 1 MeV n eq /cm 2 was used when tuning the model, the other values follow as a prediction. The error bars are calculated by varying the temperature within its estimated uncertainty combined with the uncertainty given by the deviation between the 2D and 3D simulation. to occur [22]. As the bias voltage increases, the regions with a field high enough to cause significant impact ionisation extend further into the device, and the overall field strength increases; and consequently the current caused by avalanche generation will increase. This is shown in Fig. 16(a). Fig. 16(b) shows a comparison of simulations with and without avalanche generation, confirming that the non-saturating behaviour of the I-V curve is caused by avalanche generation. Avalanche multiplication also explains why the temperature sensitivity of the charge collection efficiency increases with fluence and bias voltage, as the avalanche generation rate is temperature dependent. If charge multiplication is turned off in the simulation, the error bars due to the uncertainty in temperature are reduced by a factor of two.

Conclusions
The proposed radiation damage model for Sentaurus TCAD is able to reproduce the I-V curves of -type silicon sensors up to a fluence of 8 × 10 15 1 MeV n eq /cm 2 . The model has been tested for temperatures between −38 • C and −31 • C on a range of sensors with different irradiation types and profiles. The model captures the transition from a linear electric field and saturating I-V curve to a double junction electric field and a non-saturating I-V curve. The latter is shown to be a consequence of avalanche generation in the high-field regions of the double junction profile.
Furthermore, it is shown that a two-dimensional approximation of a pixel detector for CCE simulations is acceptable if the simulated particle passes close to the centre of the pixel. The CCE calculated using the proposed model is in agreement (within the estimated range of uncertainty) with experimental data at fluences of 4 × 10 15 1 MeV n eq /cm 2 and 8 × 10 15 1 MeV n eq /cm 2 .
It is the hope of the authors that the proposed model contribute towards obtaining a comprehensive TCAD model of radiation damage in silicon sensors. It seems clear that avalanche multiplication becomes increasingly important at high fluences, and that it must be carefully included in order to obtain a complete model.