Modelling radiation damage to pixel sensors in the ATLAS detector

TL;DR

This paper develops a detailed digitization model to simulate radiation damage effects in ATLAS pixel sensors, validated against collision data, to improve understanding of sensor performance under high radiation.

Contribution

It introduces a novel digitization model that incorporates radiation damage effects, validated with real collision data for the first time.

Findings

Model accurately predicts charge collection efficiency.

Model reproduces Lorentz angle variations.

Validated against data up to $10^{15}$ 1 MeV ${n}_{eq}/{cm}^2$.

Abstract

Silicon pixel detectors are at the core of the current and planned upgrade of the ATLAS experiment at the LHC. Given their close proximity to the interaction point, these detectors will be exposed to an unprecedented amount of radiation over their lifetime. The current pixel detector will receive damage from non-ionizing radiation in excess of $10^{15}$ 1 MeV ${n}_{eq}/{cm}^2$, while the pixel detector designed for the high-luminosity LHC must cope with an order of magnitude larger fluence. This paper presents a digitization model incorporating effects of radiation damage to the pixel sensors. The model is described in detail and predictions for the charge collection efficiency and Lorentz angle are compared with collision data collected between 2015 and 2017 ($\leq 10^{15}$ 1 MeV ${n}_{eq}/{cm}^2$).

Figures

Tables

Table 1. Table 1: Introduction rates of the Hamburg model as obtained by adjusting the simulated depletion voltage to the available measurements. For comparison, in the last column the values reported by the ROSE Collaboration  [ 6 ] are listed for oxygen-enriched silicon, separately for protons ( p 𝑝 p ) and neutrons ( n 𝑛 n ).

Parameter	IBL [$\times {10}^{-2}{\text{cm}}^{-1}$]	$B$-layer [$\times {10}^{-2}{\text{cm}}^{-1}$]	ROSE Coll. [$\times {10}^{-2}{\text{cm}}^{-1}$]
$g_{\text{A}}$	$1.4\pm 0.5$	$1.4\pm 0.5$	1.4 ($n$)
$g_{\text{Y}}$	$6.0\pm 1.6$	$6.0\pm 1.6$	2.3 ($p$), 4.8 ($n$)
$g_{\text{C}}$	$1.1\pm 0.3$	$0.45\pm 0.1$	0.53 ($p$), 2.0 ($n$)

Table 2. Table 2: Nominal predictions from the Hamburg model for the effective doping concentration N eff subscript 𝑁 eff N_{\text{eff}} and for donor (acceptor) concentration N D(A) subscript 𝑁 D(A) N_{\text{D(A)}} for two points in time during Run 2. The value of N D subscript 𝑁 D N_{\text{D}} was chosen to be numerically small (for technical reasons, it cannot be exactly zero) and the actual value has little impact on the result. The fluence 2 × 10 14 absent superscript 10 14 \times 10^{14} n / eq cm 2 {}_{\text{eq}}/\text{cm}^{2} was reached near a time of annealing where the effective doping concentration changed by about 4% over a short period in fluence. The reported doping concentration and corresponding bias voltage correspond to approximately the midpoint of the concentration during this brief period.

$\Phi$ [n${}_{\text{eq}}{}^{}/{\text{cm}}^2$]	Approx. date	$N_{\text{eff}}$ [cm-3]	$N_{\text{eff}}$ uncert. [%]	$N_{\text{D}}$ [cm-3]	$N_{\text{A}}$ [cm-3]	$V_{\text{depl}}$ [V]
1$\times {10}^{14}$	9/7/2016	$-1.62\times {10}^{12}$	9	0.02$\times {10}^{12}$	$-1.64\times {10}^{12}$	50
2$\times {10}^{14}$	8/9/2016	$-2.72\times {10}^{12}$	21	0.02$\times {10}^{12}$	$-2.74\times {10}^{12}$	85

Table 3. Table 3: Values used in TCAD simulations for deep acceptor (donor) defect concentrations N A subscript 𝑁 A N_{\text{A}} ( N D subscript 𝑁 D N_{\text{D}} ) and for their electron (hole) capture cross sections ( σ e , h A,D superscript subscript 𝜎 𝑒 ℎ A,D \sigma_{e,h}^{\text{A,D}} ) for three different fluences. Values are derived from the Chiochia model  [ 44 ] for temperature T = − 10 ∘ 𝑇 superscript 10 T=-10^{\circ} C. Reference  [ 44 ] gives values for Φ = 0.5 × 10 14 Φ 0.5 superscript 10 14 \Phi=0.5\times 10^{14} , 2 × 10 14 2 superscript 10 14 2\times 10^{14} , and 5.9 × 10 14 ​ n eq / cm 2 5.9 superscript 10 14 subscript n eq superscript cm 2 5.9\times 10^{14}\leavevmode\nobreak\ \text{n}_{\text{eq}}/\text{cm}^{2} . In between the reported values, the interpolated value is given by the average of the neighbouring low and high fluence points scaled to the target fluence: N A/D ​ ( Φ ) = 1 2 ​ ( N A/D ​ ( Φ low ) / Φ low + N A/D ​ ( Φ high ) / Φ high ) ​ Φ ≡ g int ​ Φ subscript 𝑁 A/D Φ 1 2 subscript 𝑁 A/D subscript Φ low subscript Φ low subscript 𝑁 A/D subscript Φ high subscript Φ high Φ subscript 𝑔 int Φ N_{\text{A/D}}(\Phi)=\frac{1}{2}(N_{\text{A/D}}(\Phi_{\text{low}})/\Phi_{\text{low}}+N_{\text{A/D}}(\Phi_{\text{high}})/\Phi_{\text{high}})\Phi\equiv g_{\text{int}}\Phi , where g int subscript 𝑔 int g_{\text{int}} is the effective introduction rate. For fluences below 0.5 × 10 14 0.5 superscript 10 14 0.5\times 10^{14} or above 5 × 10 14 5 superscript 10 14 5\times 10^{14}   n eq / cm 2 subscript n eq superscript cm 2 \text{n}_{\text{eq}}/\text{cm}^{2} , the value is scaled, based on the nearest reported value: g int = N A/D ​ ( Φ bench ) / Φ bench subscript 𝑔 int subscript 𝑁 A/D subscript Φ bench subscript Φ bench g_{\text{int}}=N_{\text{A/D}}(\Phi_{\text{bench}})/\Phi_{\text{bench}} , where Φ bench subscript Φ bench \Phi_{\text{bench}} is the nearest reported fluence.

$\Phi$ $\left[{\text{n}}_{\text{eq}}/{\text{cm}}^2\right]$	$N_{\text{A}}\times {10}^{-15}$ $\left[{\text{cm}}^{-3}\right]$	$N_{\text{D}}\times {10}^{-15}$ $\left[{\text{cm}}^{-3}\right]$	${\sigma }_e^{\text{A/D}}\times {10}^{15}$ $\left[{\text{cm}}^2\right]$	${\sigma }_h^{\text{A}}\times {10}^{15}$ $\left[{\text{cm}}^2\right]$	${\sigma }_h^{\text{D}}\times {10}^{15}$ $\left[{\text{cm}}^2\right]$	$g_{\text{int}}^{\text{A}}$ $\left[{\text{cm}}^{-1}\right]$	$g_{\text{int}}^{\text{D}}$ $\left[{\text{cm}}^{-1}\right]$
$1\times {10}^{14}$	$0.36$	$0.5$	$6.60$	$1.65$	$6.60$	3.6	5
$2\times {10}^{14}$	$0.68$	$1$	$6.60$	$1.65$	$6.60$	3.4	5
$5\times {10}^{14}$	$1.4$	$3.4$	$6.60$	$1.65$	$1.65$	2.8	6.8

Table 4. Table 4: Results for average and RMS of the space-charge density over the sensor bulk from TCAD simulation for different scenarios at fluences Φ = 1 × 10 14 Φ 1 superscript 10 14 \Phi=1\times 10^{14} and 2 × 10 14 2 superscript 10 14 2\times 10^{14}  n / eq cm 2 {}_{\text{eq}}/\text{cm}^{2} ; the last row corresponds to the TCAD with eff. annealing scenario. Bias voltage V bias subscript 𝑉 bias V_{\text{bias}} was 150 V. Since ρ 𝜌 \rho is constant for the Hamburg scenario, the RMS is zero.

$\Phi =1\times {10}^{14}$ n${}_{\text{eq}}{}^{}/{\text{cm}}^2$
Scenario	$⟨\rho /e⟩$ [cm-3]	RMS [cm-3]
Chiochia	$-1.5\times {10}^{12}$	0.7$\times {10}^{12}$
Hamburg	$-1.62\times {10}^{12}$	0
$N_{\text{A}}^{\text{TCAD}}$+3%	$-1.7\times {10}^{12}$	0.8$\times {10}^{12}$

Table 5. Table 5: Measurements of the trapping constant β 𝛽 \beta are summarized, normalized to a temperature of 0 ∘ C. Some measurements are reported after annealing to the minimum in the full depletion voltage V depl subscript 𝑉 depl V_{\text{depl}} (reached in about 80 minutes at 60 ∘ C) while others correspond to the asymptotic values observed after long annealing times.

Irradiation	Annealing	${\beta }_e$ (${10}^{-16}{\text{cm}}^2$/ns)	${\beta }_h$ (${10}^{-16}{\text{cm}}^2$/ns)	Reference	Method
Neutrons	minimum $V_{\text{depl}}$	$4.0\pm 0.1$	$5.7\pm 0.2$	[49]	TCT
Pions	minimum $V_{\text{depl}}$	$5.5\pm 0.2$	$7.3\pm 0.2$	[49]	TCT
Protons	minimum $V_{\text{depl}}$	$5.13\pm 0.16$	$5.04\pm 0.18$	[50]	TCT
Neutrons	$>50$ hours at 60∘C	$2.6\pm 0.1$	$7.0\pm 0.2$	[49]	TCT
Protons	$>10$ hours at 60∘C	$3.2\pm 0.1$	$5.2\pm 0.3$	[50]	TCT
Protons	minimum $V_{\text{depl}}$	$4.0\pm 1.4$	-	[3, 51]	Test-beam
Protons	25h at 60∘C	$2.2\pm 0.4$	-	[3, 51]	Test-beam

Table 6. Table 6: Values used in TCAD simulations for deep defect energies, introduction rates ( g int subscript 𝑔 int g_{\text{int}} ) and for their electron and hole capture cross sections ( σ e , h subscript 𝜎 𝑒 ℎ \sigma_{e,h} ). Values are taken from the Perugia model  [ 53 ] .

Type	Energy [eV]	$g_{\text{int}}\left[{\text{cm}}^{-1}\right]$	${\sigma }_e\left[{\text{cm}}^2\right]$	${\sigma }_h\left[{\text{cm}}^2\right]$
Acceptor	$E_{\text{c}}-0.42$	1.613	$1\times {10}^{-15}$	$1\times {10}^{-14}$
Acceptor	$E_{\text{c}}-0.46$	0.9	$3\times {10}^{-15}$	$3\times {10}^{-14}$
Donor	$E_{\text{v}}+0.36$	0.9	$3.23\times {10}^{-13}$	$3.23\times {10}^{-14}$

Table 7. Table 7: List of systematic uncertainties considered in the simulation and their relative impact on the predicted charge collection efficiency. Blank spaces correspond to uncertainties that are below 0.005 0.005 0.005 . The total uncertainty is the sum in quadrature of the maximum of the up and down variations.

Bias voltage [V]	$80$	$80$	$80$	$150$	$150$	$150$
Fluence $\left[{10}^{14}{\text{n}}_{\text{eq}}/{\text{cm}}^2\right]$	$1$	$2$	$5$	$1$	$2$	$5$
Variation	Impact [%]	Impact [%]	Impact [%]	Impact [%]	Impact [%]	Impact [%]
Energy acceptor $+10\%$	$0.4$	$2.2$	$0.7$	$0.2$	$1.6$	$1.5$
Energy donor $+10\%$	$0.5$		$4.5$	$0.2$	$0.3$	$5.7$
Energy acceptor $-10\%$		$1.7$	$3.8$	$0.1$	$0.3$	$1.6$
Energy donor $-10\%$	$0.5$	$0.1$	$2.9$	$0.1$	$0.9$
$g_{\text{int}}$ acceptor $+10\%$	$0.4$	$0.3$	$3.8$	$0.1$	$0.9$	$1.3$
$g_{\text{int}}$ donor $+10\%$	$0.2$	$1.0$	$4.2$	$0.1$	$0.4$	$5.7$
$g_{\text{int}}$ acceptor $-10\%$	$0.3$	$1.7$	$14.0$	$0.1$	$0.3$	$1.5$
$g_{\text{int}}$ donor $-10\%$	$0.03$	$0.4$	$6.8$	$0.1$	$0.7$	$6.9$
${\sigma }_e$ acceptor $+10\%$	$0.3$	$1.4$	$0.9$	$0.06$	$0.4$	$2.3$
${\sigma }_e$ donor $+10\%$	$0.2$	$1.0$	$1.8$	$0.01$	$0.4$	$0.6$
${\sigma }_e$ acceptor $-10\%$	$0.4$	$0.3$	$1.9$	$0.1$	$0.8$	$4.7$
${\sigma }_e$ donor $-10\%$	$0.1$	$0.6$	$0.1$	$0.1$	$0.6$	$5.5$
${\sigma }_h$ acceptor $+10\%$	$0.3$	$1.5$	$1.3$	$0.1$	$0.8$	$2.6$
${\sigma }_h$ donor $+10\%$	$0.3$	$0.1$	$1.4$	$0.1$	$0.6$	$5.0$
${\sigma }_h$ acceptor $-10\%$	$0.3$	$2.0$	$1.0$	$0.1$	$0.3$	$2.2$
${\sigma }_h$ donor $-10\%$	$0.3$	$0.5$	$0.8$	$0.01$	$0.3$	$0.8$
electron trapping constant $+33\%$	$2.3$	$1.1$	$1.1$	$1.3$	$2.5$	$1.0$
hole trapping constant $+23\%$	$0.5$	$2.4$	$11.8$	$1.2$	$2.1$	$6.0$
electron trapping constant $-33\%$	$1.1$	$2.0$	$5.1$	$0.5$	$0.1$	$6.1$
hole trapping constant $-23\%$	$2.2$	$2.5$	$0.4$	$0.6$	$0.2$	$2.9$
Total Uncertainty	$3.4$	$5.1$	$21.3$	$1.9$	$4.2$	$15.6$

Equations

$$I_{\text{leak}}=(\Phi /L_{\text{int}})\cdot \sum _{i=1}^n{V}_i\cdot L_{\text{int,i}}\cdot \left[{\alpha }_I\exp \left(-\sum _{j=i}^n\frac{t_j}{\tau (T_j)}\right)+{\alpha }_0^{\ast }-\beta \log \left(\sum _{j=i}^n\frac{\Theta (T_j)\cdot t_j}{t_0}\right)\right],$$

$$\Theta (T)=\exp \left[-\frac{E_I^{\ast }}{k_{\text{B}}}\left(\frac{1}{T}-\frac{1}{T_{\text{ref}}}\right)\right],$$

$$N_{\text{eff}}(t)=N_{\text{D}}^{\text{non-removable}}(0)+N_{\text{D}}^{\text{removable}}(t)-N_{\text{A}}^{\text{stable}}(t)-N_{\text{A}}^{\text{beneficial}}(t)-N_{\text{A}}^{\text{reverse}}(t),$$

$$\frac{\text{d}}{\text{d}t}{N}_{\text{D}}^{\text{removable}}(t)$$

$$\frac{\text{d}}{\text{d}t}{N}_{\text{A}}^{\text{stable}}(t)$$

$$\frac{\text{d}}{\text{d}t}{N}_{\text{A}}^{\text{beneficial}}(t)$$

$$\frac{\text{d}}{\text{d}t}{N}_{\text{N}}^{\text{reverse}}(t)$$

$$\frac{\text{d}}{\text{d}t}{N}_{\text{A}}^{\text{reverse}}(t)$$

$$V_{\text{depl}}=|N_{\text{eff}}|\cdot \frac{e{d}^2}{2ϵ{ϵ}_0},$$

$$⟨\rho /e{⟩}_{\text{TCAD}}=(N_{\text{eff}}{)}_{\text{Hamburg}}.$$

$$t_{\text{collection}}({\vec{x}}_{\text{initial}})={\int }_C\frac{\text{d}s}{r\mu (E)E},$$

$${\vec{x}}_{\text{trap}}(t_{\text{to trap}})={\int }_{\text{0}}^{\text{tto trap}}r\mu (E)\vec{E}\text{d}t,$$

$$\tan {\theta }_{\text{L}}^{\text{integrated}}(z_{\text{initial}},z_{\text{final}})=\frac{rB}{|z_{\text{final}}-z_{\text{initial}}|}{\int }_{z_{\text{initial}}}^{z_{\text{final}}}\mu (E(z))\text{d}z.$$

$$i=q{\vec{E}}_{\text{w}}\cdot \vec{v},$$

$$Q_{\text{induced}}=-q[{\varphi }_{\text{w}}({\vec{x}}_{\text{end}})-{\varphi }_{\text{w}}({\vec{x}}_{\text{start}})],$$

$$F(\alpha )=[a\times |\tan \alpha -\tan {\theta }_{\text{L}}|+b/\sqrt{\cos \alpha }]\otimes G(\alpha |\mu =0,\sigma ),$$

Full text

\AtlasTitle

Modelling radiation damage to

pixel sensors in the ATLAS detector

\AtlasJournalRefJINST 14 (2019) P06012 \AtlasDOI10.1088/1748-0221/14/06/P06012 \PreprintIdNumberCERN-EP-2019-061 \AtlasAbstractSilicon pixel detectors are at the core of the current and planned upgrade of the ATLAS experiment at the LHC. Given their close proximity to the interaction point, these detectors will be exposed to an unprecedented amount of radiation over their lifetime. The current pixel detector will receive damage from non-ionizing radiation in excess of $10^{15}$ 1 MeV $\text{n}_{\text{eq}}/\text{cm}^{2}$, while the pixel detector designed for the high-luminosity LHC must cope with an order of magnitude larger fluence. This paper presents a digitization model incorporating effects of radiation damage to the pixel sensors. The model is described in detail and predictions for the charge collection efficiency and Lorentz angle are compared with collision data collected between 2015 and 2017 ($\leq 10^{15}$ 1 MeV $\text{n}_{\text{eq}}/\text{cm}^{2}$).

\size@chapter\sectfont

Contents

@afterheading@starttoc

toc

1 Introduction

As the subdetector in closest proximity to the interaction point, the ATLAS pixel detector will be exposed to an unprecedented amount of radiation over its lifetime. The modules comprising the detector are designed to be radiation tolerant, but their performance will still degrade over time. It is therefore of crucial importance to model the impact of radiation damage for an accurate simulation of charged-particle interactions with the detector and the reconstruction of their trajectories (tracks). Modelling radiation damage effects is especially relevant for the high-luminosity (HL) upgrade of the Large Hadron Collider (LHC); the instantaneous and integrated luminosities will exceed current values by factors of 5 and 10, respectively. The simulations for the present (Run 1: 2010-12, Run 2: 2015-18) and future ATLAS detectors currently do not model the effect of silicon sensor radiation damage [1, 2].

This article documents the physics and validation of the pixel radiation damage models that will be incorporated into the ATLAS simulation. Section 2 briefly introduces the specifications of the ATLAS pixel detector and provides an overview of the impact of radiation damage effects. Measurements of the fluence and depletion voltage are presented in Section 3. A model of charge deposition and measurement that includes radiation damage effects is documented in Section 4. Comparisons and validation of the simulation with data are presented in Section 5 and conclusions are given in Section 6.

2 The ATLAS pixel detector and radiation damage effects

The ATLAS pixel detector [3, 4, 5] consists of four barrel layers and a total of six disc layers, three at each end of the barrel region. The four barrel layers are composed of $n^{+}$-in-$n$ planar oxygenated [6, 7] silicon sensors at radii of 33.5, 50.5, 88.5, and 122.5 mm from the geometric centre of the ATLAS detector [8]. The sensors on the innermost barrel layer (the insertable $B$-layer or IBL [4, 5], installed between Runs 1 and 2) are $200$ µm thick, while the sensors in the other layers are $250$ µm thick. At high $|z|$111ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the $z$-axis coinciding with the axis of the beam pipe. The $x$-axis points from the IP towards the centre of the LHC ring, and the $y$-axis points upward. Cylindrical coordinates ($r$,$\phi$) are used in the transverse plane, $\phi$ being the azimuthal angle around the $z$-axis. The pseudorapidity is defined in terms of the polar angle $\theta$ as $\eta=-\ln\tan(\theta/2)$. on the innermost barrel layer, there are $n^{+}$-in-$p$ 3D sensors [9] that are $230$ µm thick. The innermost barrel layer pixel pitch is $50\times 250$ µm2; everywhere else the pixel pitch is $50\times 400$ µm2. Charged particles traversing the sensors deposit energy by ionizing the silicon bulk; for typical LHC energies, such particles are nearly minimum-ionizing particles (MIP). The deposited charge drifts through the sensor and the analogue signal recorded by the electrode is digitized, buffered, and read out using an FEI4 [10] (IBL) or FEI3 [3] (all other layers) chip. Non-ionizing interactions from heavy particles and nuclei lead to radiation damage, which modifies the sensor bulk and can therefore alter the detection of MIPs.

Radiation damage in the sensor bulk is caused primarily by displacing a silicon atom out of its lattice site resulting in a silicon interstitial site and a leftover vacancy (Frenkel pair) [11, 12]. These primary defects build, depending on the recoil energy, cluster defects and point defects in the silicon lattice that cause energy levels in the band gap. When activated and occupied, these states lead to a change in the effective doping concentration, a reduced signal collection efficiency due to charge trapping, and an increase in the sensor leakage current that is proportional to the fluence received. The change in effective doping concentration has consequences for the depletion voltage and electric field profile. For the pixel planar sensors before irradiation, the depletion region grows from the back side of the sensor towards the pixel $n^{+}$ implant. After irradiation, the effective doping concentration decreases with increasing fluence until the sensor bulk undergoes space-charge sign inversion (often called type inversion) from $n$-type to $p$-type. After this type inversion, the depletion region grows from the $n^{+}$ implant towards the back side of the sensor and the depletion voltage gradually increases with further irradiation (more details are given in Section 3.2). The effective doping concentration is further complicated by annealing in which new defects are formed or existing defects dissociate due to their thermal motion in the silicon lattice [11]. Consequently, radiation damage effects depend on both the irradiation and temperature history. The silicon bulk of the IBL planar sensors underwent type inversion after about 3 fb*-1* of data collected in 2015 and the second innermost layer ($B$-layer) inverted in the 2012 run after about 5 fb*-1*. The outer two layers inverted between Runs 1 and 2.

In ATLAS, complex radiation fields are simulated by propagating inelastic proton–proton interactions, generated by Pythia 8 [13, 14] using the MSTW2008LO parton distribution functions [15] and the A2 set of tuned parameters [16], through the ATLAS detector material using the particle transport code FLUKA [17, 18]. Particles are transported down to an energy of 100 keV, except for photons (30 keV) and neutrons (thermal). It is important to model as accurately as possible all the inner detector and calorimeter geometry details because high-energy hadron cascades in the material lead to increased particle fluences in the inner detector, especially neutrons. A description of the ATLAS FLUKA simulation framework can be found in Ref. [19].

Predictions of the 1 MeV neutron-equivalent fluences222For silicon sensors the relevant measure of the radiation damage is the non-ionizing energy loss (NIEL), normally expressed as the equivalent damage of a fluence of 1 MeV neutrons ($\text{n}_{\text{eq}}/\text{cm}^{2}$). per fb*-1* for silicon in the ATLAS FLUKA inner detector geometry are shown in Figure 1(a). The dominant contribution is from charged pions originating directly from the proton–proton collisions. The fluence values averaged over all barrel modules for the four pixel layers starting from the innermost one are $6.1\times 10^{12}$, $2.9\times 10^{12}$, $1.2\times 10^{12}$ and $7.8\times 10^{11}$ $\text{n}_{\text{eq}}/\text{cm}^{2}/\text{fb}^{-1}$, respectively. The fluence depends on the $z$ position as the material and particle composition are $\eta$-dependent. For example, in the IBL the maximum predicted value of $6.6\times 10^{12}$ $\text{n}_{\text{eq}}/\text{cm}^{2}/\text{fb}^{-1}$ in the central location is about $10\%$ higher than in the end regions (studied further in Section 3.1). Figure 1(b) shows the 1 MeV neutron-equivalent fluence as a function of time, based on the FLUKA simulation. The luminosity is determined by a set of dedicated luminosity detectors [20] that are calibrated using the van der Meer beam-separation method [21]. By the end of the proton–proton collision runs in 2017, the IBL and $B$-layer had received integrated fluences of approximately $\Phi=6\times 10^{14}$ and $3\times 10^{14}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$, respectively. The two outer layers have been exposed to less than half the fluence of the inner layers.

The goal of this paper is to present a model for radiation damage to silicon sensors that is fast enough to be incorporated directly into the digitization step of the ATLAS Monte Carlo (MC) simulation, i.e. the conversion from energy depositions from charged particles to digital signals sent from module front ends to the detector read-out system. In the context of the full ATLAS simulation chain [1], digitization occurs after the generation of outgoing particles from the hard-scatter collision and the simulation of their interactions with the detector and before event reconstruction, which is the same for data and simulation. The CMS Collaboration has developed a model of radiation damage [22, 23, 24, 25, 22],333This model is used in some HL-LHC projection studies, but there is currently no public documentation with a detailed description of the implementation in the CMS software. validated with test-beam data, but it is used to apply template corrections to the total deposited charge in simulation from a model without inherent radiation damage effects and so is not directly comparable to the methods described here.

There are two types of microscopically motivated effective radiation damage models used for the studies presented here: Hamburg444See Ref. [11] and references therein. This model is a phenomenological approach that includes some physically well-motivated components and other aspects that are not directly based on the microphysics of defect states. and models developed in the framework of Technology Computer Aided Design (TCAD) simulations. In reality, the microphysics is complex, involving many defect states, but each model includes a small number of effective components to capture the main effects. The Hamburg model describes annealing and is only used to validate conditions data (Section 3). Stand-alone implementations of this model simulate the time-dependent leakage current (Section 3.1) and doping concentration (Section 3.2) for checking the fluence and depletion voltage. The second type of model (TCAD) is used directly in the digitizer (software that performs digitization) described in Section 4. In contrast to the Hamburg model, radiation damage implemented in TCAD predicts a non-uniform spatial distribution of space-charge density and thus a more realistic electric field profile (Section 4.2) for computing charge propagation inside the sensor bulk.

Multiple radiation damage models are required since no model accounts for both annealing and a non-uniform space-charge density distribution. Therefore, each model is used where it is most appropriate. An approximate combination of model predictions is described in Section 4.2.4. However, for the present levels of annealing, the combination yields variations in electric field profiles that are smaller than the uncertainty in the TCAD radiation damage model parameters (Section 4.2.3) and so is not used for the final results (Section 5) – only TCAD input without annealing is currently used for the digitizer. In the future, when there is more annealing and the radiation damage model parameters are further constrained from data, it will become a crucial and challenging project to combine the power of both types of models.

3 Validating sensor conditions

3.1 Luminosity to fluence

The most important input to the radiation damage digitization model is the estimated fluence. Section 2 introduced the baseline FLUKA simulation that is used to determine the conversion factor ($\Phi/L_{\text{int}}$) between integrated luminosity and fluence. In order to estimate systematic uncertainties in these predictions, the fluence is converted into a prediction for the leakage current. The leakage current can be precisely measured and therefore provides a solid validation for the FLUKA simulation. For $n$ time intervals, the predicted leakage current is given by Ref. [11]:

[TABLE]

where L_{\text{int,i}} is the integrated luminosity, $t_{i}$ is the duration, and $T_{i}$ is the temperature in time interval $i$. The first sum is over all time periods and the two sums inside the exponential and logarithm functions are over the time between the irradiation in time period $i$ and the present time. The other symbols in Eq. (1) are $t_{0}=1\;$min, $V_{i}$ is the depleted volume (in cm3), $\alpha_{I}=(1.23\pm 0.06)\times 10^{-17}\text{ A}/\text{cm}$, $\tau$ follows an Arrhenius equation $\tau^{-1}=(1.2^{+5.3}_{-1.0})\times 10^{13}\text{ s}^{-1}\times\text{e}^{(-1.11\pm 0.05)\text{ eV}/k_{\text{B}}T}$, where $k_{\text{B}}$ is the Boltzmann constant, $\alpha_{0}^{*}=7.07\cdot 10^{-17}\;$A/cm, and555A small temperature dependence has been observed in the value of $\beta$ [11]. For this analysis, the reported value at 21∘C is taken as it is closest to the operational temperature range of the detector. $\beta=(3.29\pm 0.18)\times 10^{-18}$ A/cm. The time scaling function $\Theta(T)$ is defined by666This is not the only way to incorporate time-dependence in the thermal history. Another proposal is to sum the inverse temperatures [26]. Such a method has been compared with Eq. (2) and results in similar predictions for the leakage current with the current fluence levels and annealing times.

[TABLE]

where $E_{I}^{*}=(1.30\pm 0.14)$ eV and $T_{\text{ref}}$ is a reference temperature, typically 20∘C.

Using the measured module temperature as a function of time, Eq. (1) is used to predict the leakage current as shown in Figure 2. The leakage current is scaled to correspond to a temperature of 20∘C using the factor (see e.g., Ref. [27]) $(T_{\text{R}}/T)^{2}\exp\left(-E_{\text{eff}}(T_{\text{R}}^{-1}-T^{-1})/2k_{\text{B}}\right)$, where $T_{\text{R}}=20^{\circ}$C and $E_{\text{eff}}=1.12$ eV. The value of $E_{\text{eff}}$ is lower than the one measured in Ref. [28], but was found to agree better with the data. Measurements of the properties describing the modules were updated every ten minutes. Since the IBL was newly inserted before the 2015 run, the initial leakage current level is compatible with zero. A constant $\Phi/L_{\text{int}}$ conversion factor is fit to the data per module group. Module groups differ by their distance along the beam direction from the geometric centre of the detector. Each module group is $8$ cm long on both sides of the detector along the beam direction. The groups M1, M2, M3 approximately span the ranges $z\in[-8,8]$ cm, $|z|\in[8,16]$ cm, and $|z|\in[16,24]$ cm, respectively. The M4 modules use 3D sensors; M4 spans the range $|z|\in[24,32]$ cm. Only the time interval indicated by a dashed region in Figure 2 is used in the fluence rate extraction. Prior to this time, the IBL was under-depleted and after this time, the frequency of measurements decreased. A sensor volume correction is applied to the under-depleted data. After this correction, the adjusted simulation reproduces the trends observed in the data both inside and outside of the fit region.

Figure 3 shows that there is a stronger $z$-dependence in the measured fluence compared with the FLUKA predictions described in Section 2. The error bars on $\Phi/L_{\text{int}}$ predicted by the Hamburg model fitted to data are dominated by a conservative $10\%$ uncertainty, accounting for the possible difference between the leakage current at the operational bias voltage and the current at the full depletion voltage (see Section 3.2). After irradiation, the leakage current increases with increasing bias voltage also after full depletion, while the Hamburg model predicts a constant leakage current above the full depletion voltage. Therefore, the choice of voltage for the leakage current measurement is crucial for comparison with the Hamburg model prediction [29]. Uncertainties due to the annealing model ($0.1\%$) and data fit ($0.5\%$) are subdominant. The predictions in Figure 3 deviate from the measured values by about $1.5\sigma$ of the uncertainty at $z=0$, with larger deviations at higher $|z|$. In addition to the Pythia+FLUKA prediction described in Section 2, Figure 3 also shows predictions with an updated Pythia set of tuned parameters (A3 [30]) as well as an alternative geometry and transport model using Geant4 [31]. Neither of these variations can account for the $z$-dependence, but this does illustrate part of the uncertainty due to the transport model and particle generator. There is also a significant source of uncertainty from the silicon hardness factors [12] (common to both the Geant4 and FLUKA models). The hardness factors used here are from the RD50 database [32, 33, 34, 35, 36], but all of these values are without uncertainty and many are based only on simulation. The uncertainty in the hardness factors affects both the prediction and the Hamburg model (through the $\alpha$ parameters). Future collision data may be able to constrain these hardness factors. As shown in Fig. 3, most of the damage is due to charged pions, protons, and neutrons, so the larger uncertainties on other particle species is a subdominant source of total uncertainty for the hardness factors.

The remainder of this paper focuses on central $|z|\approx 0$, using the FLUKA simulations without modification for the central value, but with a $15\%$ uncertainty in the fluence taken from this leakage current study. The ATLAS tracking acceptance is $|\eta|<2.5$, which corresponds to $|z|<20$ cm in the IBL.

3.2 Annealing and depletion voltage

As already introduced in Section 2, the irradiation and thermal history are accounted for in the prediction of the effective doping concentration with the Hamburg model. In this model, the effective doping concentration ($N_{\text{eff}}(t)$) has the following form:

[TABLE]

where $N_{\text{D}}^{\text{(non)-removable}}(0)$ is the initial concentration of (non)-removable donors777Where $N_{\text{eff}}(0)=N_{\text{D}}^{\text{non-removable}}(0)+N_{\text{D}}^{\text{removable}}(0)$. and the other terms are described below. The fraction of removable donors at the doping concentrations used for silicon sensors is predicted to be 100% of the initial doping concentration for charged-particle irradiation, which dominates the inner pixel layers in the ATLAS detector. The time-dependence of the terms on the right-hand side of Eq. (3) are described by the following differential equations:

[TABLE]

where $\phi(t)$ is the irradiation rate in $\text{n}_{\text{eq}}/\text{cm}^{2}/\text{s}$. Equation (4) represents the effective removal of the initial donors by mobile defects. The removal constant is $c=6.4\times 10^{-14}\text{ cm}^{2}$ [11]. The second equation, Eq. (5), represents the constant addition of stable (non-annealable) defects which act electrically as acceptors. Two additional defects are introduced in Eqs. (6) and (7). These defects, introduced during irradiation with introduction rates $g_{\text{A}}$ and $g_{\text{Y}}$, are short-lived at sufficiently high temperatures ($\gtrsim 10^{\circ}$C). The temperature-dependence of the decay rates is modelled with an Arrhenius equation, $k_{i}(T)=k_{i,0}\,\text{e}^{-E_{i}/k_{\text{B}}T}$, where $k_{\text{A,0}}=2.4^{+1.2}_{-0.8}\times 10^{13}/\text{s},k_{\text{Y,0}}=1.5^{+3.4}_{-1.1}\times 10^{15}/\text{s}$, $E_{\text{A}}=(1.09\pm 0.03)$ eV and $E_{\text{Y}}=(1.33\pm 0.3)$ eV [11]. For the beneficial annealing (Eq. (6)), the acceptor-like defects introduced during irradiation decay into neutral states with a time constant that is $\mathcal{O}(\text{days})$ at $20^{\circ}$C and $\mathcal{O}(\text{years})$ at $-15^{\circ}$C . In contrast, for reverse annealing, neutral defects are introduced during irradiation (Eq. (7)). The neutral defects can decay into acceptor-like states (Eq. (8)), decreasing (increasing) the effective doping concentration before (after) space-charge sign inversion. The timescale for reverse annealing is $\mathcal{O}(\text{weeks})$ at $20^{\circ}$C.

While the introduction rates $g_{\text{C}},g_{\text{A}}$, and $g_{\text{Y}}$ have been measured elsewhere (e.g. Ref. [6]), the reported values vary significantly amongst different materials and irradiation types, and so are fit with depletion voltage data from the ATLAS pixel detector. The notion of full depletion is not well-defined for highly irradiated sensors where the regions inside the sensor bulk can have a very low field (see Section 4.2.2). However, at moderate fluences, the depletion region is well-defined and is important for calibrating the parameters of the Hamburg model specifically for the ATLAS pixel sensors, as the full depletion voltage ($V_{\text{depl}}$) is calculated in terms of the effective doping concentration:

[TABLE]

where $d$ is the sensor thickness, $e$ is the charge of the electron, $\epsilon$ is the dielectric constant, and $\epsilon_{0}$ is the vacuum permittivity. Figure 4 shows the calculated $V_{\text{depl}}$ using the Hamburg model as a function of time for (a) the IBL and (b) the $B$-layer. In situ measurements of $V_{\text{depl}}$ of the sensors were performed with two different methods using the ATLAS pixel detector: the cross-talk scan and the bias voltage scan. The first method uses the cross-talk between adjacent pixels (square points in Figure 4). Since the pixels are isolated (i.e., no cross-talk) only after full depletion, this is a powerful measurement tool. However, this method is only applicable before space-charge sign inversion since afterwards the pixels are already isolated at low bias voltages, much before full depletion. The bias voltage scan uses the mean time over threshold (ToT) [37] of clusters of hits on reconstructed particle trajectories, measured in units of bunch crossings (25 ns). The depletion voltage is extracted by fitting two linear functions to the rising and plateau regions of the measured data. The intersection of the two lines is defined to be the depletion voltage (circular points in Figure 4). The initial calculated $V_{\text{depl}}$ is chosen to match the value measured during quality assurance of the IBL planar pixel sensors. The total uncertainty band for the calculations is due to varying the input parameters within their uncertainties in addition to a 20% uncertainty in the initial doping concentration (see e.g. Ref. [38]).

Due to the huge parameter space given by the defect introduction rates, the time necessary for the simulation, and a focus on physically rather than mathematically correct parameter combinations, the adjustment of the introduction rates888The initial effective doping concentration also has some uncertainty, but the fitted parameters are mostly set by the measurements following space-charge sign inversion and therefore are largely insensitive to this uncertainty. was performed using particular periods of time or available data, described in the following. The derived introduction rates are summarized in Table 1. The central value of $g_{\text{A}}$ is from the literature [6] since the measurements reported here were performed at times where there was no sensitivity to beneficial annealing. The uncertainty in $g_{\text{A}}$ reported in Table 1 is determined by adjusting $g_{\text{A}}$ according to the same prescription as for $g_{\text{Y}}$ and $g_{\text{C}}$ – one parameter is varied at a time until there is a large deviation (more details provided below). The value of $g_{\text{Y}}$ was extracted from the reverse annealing during the long shutdown 1 (LS1), which was an extended period when the detector was maintained at room temperature without further irradiation. Since the IBL was installed after LS1 and has not undergone significant reverse annealing, the $g_{\text{Y}}$ value of the $B$-layer is used also for the IBL. In contrast to $g_{Y}$ and $g_{A}$, $g_{C}$ can be well-constrained during any data-taking period when constant damage is accumulated. The extracted values for $g_{\text{C}}$ are different for the IBL and the $B$-layer because the particle compositions are different (relatively more neutrons for the $B$-layer and more charged pions for the IBL). The uncertainties arise from the procedure, from the luminosity-to-fluence conversion, and from the uncertainty in the temperature in the actual sensor.

For comparison, the values obtained by the ROSE Collaboration [6] for oxygen-enriched silicon are also reported. The values for $g_{\text{C}}$ and $g_{\text{Y}}$ are within the range given by the ROSE Collaboration when neutron and proton irradiations are considered. The predictive power of the simulation would benefit from more precise measurements of $g_{\text{A}}$, $g_{\text{Y}}$ and $g_{\text{C}}$, which may be possible with future ATLAS data, but are beyond the scope of the present study.

Table 2 collects predictions from the Hamburg model for the effective doping concentration $N_{\text{eff}}$ in the IBL for two points in time based on the parameter values discussed above, corresponding to lifetime fluence values of $1\times 10^{14}$ and $2\times 10^{14}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$, respectively. The thermal history of the IBL modules was taken into account. The uncertainty column includes all contributions shown in Figure 4. For the uncertainty in the depletion voltage fit to the introduction rates, one parameter is varied at a time until there is a large deviation; for the luminosity-to-fluence conversion uncertainty, the fluence is varied by $\pm 15\%$ (see Section 3.1), and for the temperature uncertainty, the input temperature in all phases is varied by $\pm 5^{\circ}$C. All three sources are added in quadrature to determine the total uncertainty.

The operational conditions of the sensor bulk studied in this section are crucial inputs to the simulation of digitization to be presented in Section 4. Overall, the Hamburg model provides an excellent description of the shape of the leakage current dependence on time; FLUKA + Pythia 8 predict the fluence at $|\eta|\approx 0$ within 15%, but deviate much more at higher $|z|$. Even though the Hamburg model does not incorporate a non-uniform electric field, it accurately describes the depletion voltage dependence on time at the current fluence levels. This may need to be revisited in Run 3 when there will be significant distortions in the space-charge density spatial distribution due to radiation damage.

4 Digitizer model

4.1 Overview

Figure 5 presents a schematic overview and flowchart of the physics models included in the digitization model. Upon initialization, the digitizer receives global information about the detector geometry (pixel size and type, tilt angle) and conditions, including the sensor bias voltage, operating temperature, and fluence. For the calculation of individual charge deposits within the pixel sensor, the digitizer takes as input the magnitude and location of energy deposited by a charged particle, and outputs a digitized encoding of the measured charge. The input is produced by Geant4 with possible corrections for straggling in thin silicon [39]. A TCAD tool is used to model the electric field, including radiation damage effects (Section 4.2.1).

The TCAD simulations consider a limited number (2–3) of effective deep defects that capture the modification of macroscopic variables such as leakage current, operational voltage and charge collection efficiency. Ionization energy is converted into electron–hole pairs ($\sim 3.65$ eV/pair) which experience thermal diffusion and drift in electric and magnetic fields. In order to speed up the simulation, groups of $\mathcal{O}(10)$ charge carriers drift and diffuse toward the collecting electrode (electrons) or back plane (holes), with a field- and temperature-dependent mobility. Charge groupings are chosen by dividing the deposited energy into a fixed number of pieces. The number of fundamental charges per charge grouping is set to be small enough so that the overestimation of fluctuations is negligible.999Suppose a fraction $(1-p)$ of electrons are trapped while drifting toward the electrode, assuming $p$ is constant for illustration and ignoring holes. If $n$ electrons are deposited, the number of electrons that reach the electrode is $np$ on average with a variance of $p(1-p)n$. If instead $m<n$ charge groupings with $n/m$ electrons per group are propagated and have trapping probability $p$, the average number of electrons that reach the electrode is still $mp\times(n/m)=np$ but the variance is $p(1-p)n^{2}/m\geq p(1-p)n$. For each charge group, a fluence-dependent time-to-trap (Section 4.5) is randomly generated and compared with the drift time (Section 4.3). If the drift time is longer than the time-to-trap, the charge group is declared trapped, and its trapping position is calculated. Since moving charges induce a current in the collecting electrode, a signal is induced on the electrodes also from trapped charges during their drift. This induced charge also applies to neighbouring pixels, which contributes to charge sharing. The induced charge is calculated from the initial and trapped positions using a weighting (‘Ramo’) potential (Section 4.6). The total induced charge is then converted into a ToT that is used by cluster and track reconstruction tools.

The schematic diagram in Figure 5(a) shows a planar sensor, but the digitization model also applies for 3D sensors. In the simulation, the only differences between planar and 3D pixels are that different TCAD models are used (Section 4.2.1) and charge carrier propagation occurs in two dimensions (transverse to the implants) instead of one (perpendicular to the electrodes). The digitization model description (Section 4) and validation (Section 5) focus on planar sensors, in part because they constitute most of the current ATLAS pixel detector and the 3D sensors are formally outside of the tracking acceptance ($|\eta|<2.5$). Some 3D sensor simulation results are nonetheless described in Section 4.7 in order to highlight the main differences relative to planar sensors.

4.2 Electric field

The radiation-induced states in the silicon band gap affect the electric field in the pixel cells by altering the electric field distribution in the bulk.101010There are also changes at the surface, but the focus here is on the deformations of the electric field within the sensor. Since the signal formation in silicon sensors depends directly and indirectly on the electric field shape (Sections 4.3, 4.4), a careful parameterization of the field profile is required. Section 4.2.1 introduces the default two-trap TCAD model used for subsequent studies. The resulting field profiles are shown in Section 4.2.2. In Section 4.2.3, systematic uncertainties in electric field profiles determined by TCAD simulations are discussed. This section ends with the presentation of a method to incorporate annealing in Section 4.2.4.

4.2.1 Simulation details

Since the charge collection is significantly different in planar and 3D sensors due to the different electrode geometries, two different set-ups are used to implement the radiation damage in the simulation model. The simulation is set up for both sensor types, but the focus is on planar sensors as the 3D sensors are outside of the standard $|\eta|<2.5$ tracking acceptance. Validation studies in Section 5 are therefore only presented for planar sensors. In TCAD simulations, impurities are only added and not removed;111111This is because structure simulation and device simulation are two separate processes in TCAD. See e.g. Section 3.3.1 in Ref. [42]. therefore one must balance initial shallow defects with radiation-induced defects. As a result, different TCAD models of effective defect states are used for each bulk type. Since the planar sensors are $n$-type and the 3D sensors are $p$-type, different TCAD models are used for the two sensors. Details for the 3D sensor simulation can be found in Section 4.7.

Investigations of the electric field profile in the bulk of irradiated silicon sensors have shown that the electric field is no longer linear with the bulk depth after irradiation (see, for example, Refs. [43, 44]). Irradiated planar sensors with non-linear profiles are simulated using the Chiochia model [44], implemented in the Silvaco TCAD package [45, 42]. The Petasecca [46] $n$-type model was also investigated, but was found to not predict space-charge sign inversion below $\Phi=1\times 10^{14}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$ and was therefore not considered further.

The simulation is performed over an area that corresponds to a quarter of an ATLAS IBL pixel sensor cell, to take advantage of symmetry. The electric field is computed at $T=-10^{\circ}$C using an effective doping concentration of $1.6\times 10^{12}/\text{cm}^{3}$ (corresponding to about 50 V full depletion voltage for unirradiated sensors [38]) with a discretization resolution of $1$ $\text{\textmu m}^{2}$. During Run 2, the operational temperature of the pixels was adjusted multiple times. For example, the IBL temperature was set to $-4^{\circ}$C in 2015, +20∘C for the first part of 2016, +10∘C for the rest of 2016, and was $-15^{\circ}$C in 2017. The TCAD simulations are all performed at $-10^{\circ}$C since this is where the models were developed. A naive temperature variation from scaling the trap occupation probability according to $\exp(-E_{\text{t}}/k_{\text{B}}T)$ ($E_{\text{t}}$ is the trap energy) predicts variations in the leakage current that are about 20% larger than the observations.121212TCAD simulations with the Chiochia model were performed at $\Phi=10^{14}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$ at 150 V and between $-20^{\circ}$C and 20∘C. The rescaling factor between $-10^{\circ}$C and the standard 20∘C is 20% lower with the Chiochia model compared with other studies in the literature [28]. The reason is that the TCAD models only include a small number of effective states, and in reality the temperature dependence is reduced when a more complex (but computationally intractable) combination of states is present. The trap energy level $E_{\text{t}}$ varies by 10% of thermal energy $k_{\text{B}}T$ (see Section 4.2.3) are found to be consistent with naive temperature variations that bracket all Run 2 operational temperatures ($-15^{\circ}$C to $+20^{\circ}$C) and therefore provide a conservative bound on the predictions presented in Section 5. In the future, high-statistics collision data may be used to tune models in situ and avoid this complication.

The Chiochia model is a double-trap model with one acceptor and one donor trap with activation energies set to $E_{\text{c}}-0.525$ eV and $E_{\text{v}}+0.48$ eV [43] for the conduction band energy level $E_{\text{c}}$ and the valence band energy level $E_{\text{v}}$, respectively. This model was developed using CMS diffusion oxygenated float zone $n^{+}$-in-$n$ pixel module prototypes and was chosen since the bulk material type is the same as in the ATLAS IBL and pixel layers and the initial effective doping concentration is similar: 50 V depletion voltage for the ATLAS IBL and 75 V for CMS with slightly thicker sensors. Sensor annealing in Ref. [44] is different than for the operational ATLAS detector, but the partially unaccounted annealing situation is incorporated in the model variations discussed in Section 4.2.2. Table 3 documents the radiation damage model parameter values used for the planar sensor TCAD simulations. Note that the concentrations are not comparable to the ones from the model presented in Table 2 because the defects presented here are deep traps while the ones for the depletion voltage Hamburg model are shallow and thus have a higher occupation probability.

4.2.2 Electric field profiles

For planar sensors, the field is largely independent of $x$ and $y$, and perpendicular to the sensor surface. Figure 6 shows the $z$-dependence of the electric field, averaged over $x$ and $y$, for an ATLAS IBL planar sensor for various fluences and bias voltages of 80 V and 150 V. Before irradiation the field is approximately linear as a function of depth. Just after type inversion (at about $2\times 10^{13}$ $n_{\text{eq}}/\text{cm}^{2}$), the field maximum is on the opposite side of the sensor. With increasing fluence, there is a minimum in the electric field in the centre of the sensor. For a fluence of $\Phi=5\times 10^{14}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$ and a bias voltage of 80 V (for which the sensors are not fully depleted, as shown in Section 3.2), this minimum is broad and occupies nearly a third of the sensor. The rest of the section considers bias voltages of 80 V and 150 V as they were the operational voltages of IBL planar sensors in 2015 and 2016, respectively.

4.2.3 Electric field profile uncertainties

In this section, systematic uncertainties in the electric field profiles evaluated using TCAD simulations are discussed. This includes studying other radiation damage models for TCAD simulations as well as the effect of varying the Chiochia model parameter values.

Extensive model comparisons are beyond the scope of this work, but the data presented in Section 5 can be used to constrain various simulations as well as tuning parameters and to derive systematic uncertainties for predictions for higher-luminosity data. In addition to the Chiochia model for the planar sensors, the Petasecca model [46] was also briefly investigated. While the model itself is supported by test-beam data, it is found to disagree qualitatively on the fluence for type-inversion with the Chiochia model131313Around $3\times 10^{14}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$ (Petasecca) versus $5\times 10^{13}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$ (Chiochia); the IBL inverted around $2\times 10^{13}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$ and the $B$-layer inverted around $2$–$3\times 10^{13}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$ (based on the measurements presented in Figure 4). and does not reproduce the observed trend of the Lorentz angle data as described later in Section 5.3. Therefore, this alternative model was not studied in further detail.

Next, the Chiochia model parameters are varied. Each parameter (capture cross sections and introduction rate) is varied by $\pm$10% of its value except the trap energy level $E_{\text{t}}$, which is varied by $\pm$10% of the thermal energy $V_{\text{th}}=k_{\text{B}}T$. The energy of the trap $E_{\text{t}}$ is defined as the energy difference between the trap and the relevant band (conduction for the acceptor-like trap and valence for the donor-like trap). The value 10% was chosen for illustration in the absence of experimental input; ideally future models or model tunings will provide quantitative uncertainty estimates.

Figure 7 shows the electric field for variations in the acceptor trap parameters for a fluence of $10^{14}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$ and a bias voltage of 80 V. The normalization of all the curves is fixed by the bias voltage and therefore all the curves cross at a point. Variations in the capture cross sections and introduction rate ($g_{\text{int}}$) introduce a change in the peak electric field that is between 15% and 30%. Similar variations are observed when the trap concentrations are varied by $\pm 10\%$ and the energy levels are varied by $\pm 10\%$ of the thermal energy $V_{\text{th}}$, which corresponds roughly to $0.4\%$ of the energy level. The latter number is chosen as a benchmark because the occupancy probability scales exponentially with the energy as $\sim\text{e}^{-E_{\text{t}}/k_{\text{B}}T}$ [47]. For example, when the acceptor energy is moved closer to the conduction band by $0.4\%$, the electric field looks symmetric around the mid-plane; moving the acceptor even closer to the conduction band would likely result in depletion starting from the back side.

As expected, the results for donor traps (not shown) show a behaviour that is opposite to the acceptor case when concentrations are changed. All the observed changes in the electric field are consistent with expectations (see for example Ref. [47]).

4.2.4 Effective modelling of annealing effects in TCAD simulations

There is no known recipe to include the annealing effects presented in Section 3.2 in TCAD-based predictions. One challenge for incorporating annealing effects is that both the Hamburg and TCAD models are motivated by multiple effective traps [11, 43, 44] and the effective states are not in one-to-one correspondence (in particular, no cluster defects are directly reproduced by TCAD simulations). In addition to this, the relative abundance of the measured acceptor-like traps changes with annealing. Lastly, the Hamburg model does not make a prediction for the dependence of the space-charge density on depth while the TCAD model predicts a non-trivial dependence, resulting in the complicated electric field profile discussed in Section 4.2.2. The non-constant space-charge density from the TCAD model is shown in Figure 8 for an ATLAS IBL planar sensor after radiation damage. For $\Phi=1\times 10^{14}$ n${}_{\text{eq}}/\text{cm}^{2}$ the space-charge density is negative and shows an almost linear dependence on the bulk depth, whereas for higher fluences the functional form is more complicated, exhibiting sizeable regions where the space-charge density is positive, in agreement with the model first proposed in Ref. [43]. This results in the non-trivial electric field profiles shown in Figure 6.

The non-constant space-charge density, despite the simulated traps being uniformly distributed across the sensor bulk, is due to the thermally generated electrons and holes which, drifting in opposite directions and getting trapped along their trajectory, give rise to a more negative (positive) region close to the electrode collecting the electrons (holes) [43]. Deviations from linearity in the space-charge density region distribution with respect to the position in the bulk are predicted by the TCAD simulations when the voltage is fixed and the fluence gets larger, as can be seen in Figure 8. These deviations can be understood in the following way: as the depletion regions develop from both sides, for fixed voltage and larger fluences, the mid part of the sensors is not depleted. Hence the space-charge density region profile deviates from linearity there.

One way to emulate annealing effects from the Hamburg model in the TCAD simulation is to match141414They do not agree exactly because the space-charge density in TCAD is dynamically generated and not known a priori. the effective doping concentration predictions from the former,151515The physical origin of the effective doping concentration is not exactly the same for the Hamburg and TCAD models. The approach given here is a first approximation that must be expanded upon in the future when annealing effects are much more prominent. such as the ones presented in Table 2, to the average space-charge density (normalized by the electron charge) of the latter:

[TABLE]

Two different scenarios to realize the situation described in Eq. (9) are studied: the Hamburg scenario and one in which the concentration of acceptor traps in TCAD simulations $N_{\text{A}}^{\text{TCAD}}$ was changed to satisfy Eq. (9), referred to in the following as TCAD with effective annealing. For the sake of comparison a third one was added, called the Chiochia scenario, which is the default set-up described in Section 4.2.1 with no modifications to emulate annealing.

For the Hamburg scenario, the concentration of shallow donors in the structure is set to a very low value and the deep acceptor and donor concentrations are adjusted as a function of depth in order to create a constant space-charge density as predicted by the static Hamburg model everywhere in the bulk. The Hamburg scenario is qualitatively different than the Chiochia one and would predict an electric field that is linear (more below), which is in contrast to various measurements elsewhere [44]. For the TCAD with effective annealing scenario, the space-charge density can vary in the bulk and Eq. (9) was solved by varying the acceptor concentrations $N_{\text{A}}^{\text{TCAD}}$.

Figure 9 shows the space-charge density predicted by TCAD simulations in the three scenarios with a bias voltage of 150 V and at the points in the irradiation and temperature history reported in Table 2. The average space-charge density in the various scenarios is summarized in Table 4 for the two fluences shown in Figure 9. The electric field profiles corresponding to the three scenarios shown in Figure 9 are presented in Figure 10. For the Hamburg scenario (shown only for comparison), the profile is linear while this is not the case for the other scenarios, especially at the higher fluence $\Phi=2\times 10^{14}$ n${}_{\text{eq}}/\text{cm}^{2}$.

In summary, the TCAD with effective annealing scenario used an acceptor trap density $N_{A}^{\text{TCAD}}$ in the TCAD simulations that was increased by 3% at a fluence of $\Phi=1\times 10^{14}$ n${}_{\text{eq}}/\text{cm}^{2}$ and reduced by 1.6% at a fluence $\Phi=2\times 10^{14}$ n${}_{\text{eq}}/\text{cm}^{2}$ to emulate the effect of annealing predicted by the Hamburg model. These variations are currently within the model variations described in Section 4.2.3 that are used to set systematic uncertainties on the radiation damage model parameters. Therefore, no corrections or uncertainties are applied to the simulation to account for annealing for the current radiation levels. This must be revisited when the IBL has experienced significant annealing.

4.3 Time-to-electrode, position-at-trap

Numerically, propagating charges through the silicon sensor can be computationally expensive, but fortunately can be computed once per geometry and set of conditions (temperature, bias voltage, and fluence). Electrons and holes drift with a carrier-dependent mobility ($\mu$) that depends on the electric field ($E$) and temperature [48]. The drift velocity is given by $\vec{v}(E)\sim r\mu(E)\vec{E}$ ($r=\text{Hall}$ scattering factor, $r\approx 1$) and the charge collection time is estimated via

[TABLE]

where $C$ is the path from $\vec{x}_{\text{initial}}$ to $\vec{x}_{\text{final}}$ that is determined by the equations of motion $\vec{v}=r\mu(E)\vec{E}$; $\vec{x}_{\text{final}}$ depends on the type of the charge carrier. For planar sensors, the field is nearly independent of $x$ and $y$, so the time to the electrode is parameterized in $z$ and the integral in Eq. (10) is one-dimensional. Since the mobility of holes is much lower than for electrons, it takes holes much longer (factor of $\gtrsim 3$–$5$) on average to arrive at the electrode. The collection time varies with fluence, bias voltage, and distance to the electrode, but is on average $\mathcal{O}(1)$–$\mathcal{O}(10)$ ns for $\Phi\lesssim 10^{15}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$ as shown in Figure 11.

For charges that are trapped (see Section 4.5), the location of the trapped charge must be known. The position-at-trap can be calculated in a fashion similar to the time-to-electrode from Eq. (10). In particular, the location is given by

[TABLE]

where $t_{\text{to trap}}$ is either the drift time (if the charge is not trapped) or a random time set by the trapping constant (Section 4.5). Representative position-at-trap maps are shown in Figure 12 for planar sensors. If the time travelled is larger than the drift time, then the electrons reach the collecting electrode ($z=0$) and holes reach the back side ($z=200$ µm). The corresponding maps for 3D sensors are more difficult to visualize due to their higher dimensionality.

4.4 Lorentz angle

The Lorentz angle ($\theta_{\text{L}}$) is the result of balancing electric and magnetic forces, and is defined as the incidence angle that produces the smallest cluster size in the transverse direction. As $\theta_{\text{L}}$ depends on the shape of the electric field, it is affected by radiation-induced changes of the space-charge density. Changes in the electric field affect $\theta_{\text{L}}$ indirectly through the mobility as $\tan\theta_{\text{L}}=r\mu(E)B$, where $B$ is the magnetic field. Figure 13(a) demonstrates the change in the Lorentz angle along the trajectory of electrons and holes. As the mobility increases with decreasing electric field strength, the Lorentz angle is largest near the centre of the sensors when irradiated. The path-dependence of the Lorentz angle is modelled in a manner similar to the position maps from Section 4.3 (Figure 12) by averaging the Lorentz angle along the path:

[TABLE]

The drift along the $\phi$ direction is then modified as $|z_{\text{final}}-z_{\text{initial}}|\tan\theta_{\text{L}}^{\text{integrated}}(z_{\text{initial}},z_{\text{final}})$, where the direction is the same for both the electrons and holes because both the charge and velocity sign are reversed for holes relative to electrons. The integrated Lorentz angle variations are shown in Figure 13; for this fluence and bias voltage, the integrated Lorentz angle can change by as much as a factor of two, depending on the starting and ending position.

4.5 Charge trapping

In the simulation, charge carriers are declared trapped if the projected time to reach the electrode, as defined in Section 4.3, exceeds a randomly set trapping time $t$ that is exponentially distributed with mean value $1/(\beta\Phi)$ [49, 50], where $\Phi$ is the fluence and $\beta$ is the trapping constant.

The linear relation with fluence has been measured and shown to hold with very good precision up to $10^{15}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$, but the value of $\beta$ has been found to depend on the type of irradiation, the temperature, the annealing history of the device, and the type of charge carriers (electrons or holes) [49, 50]. The measurements in Refs. [49, 50] were performed with the transient current technique (TCT). The results of the measurements are reported in Table 5. Since measurements were performed at different temperatures between $-10^{\circ}$C and 10*∘C and a significant decrease of $\beta$ with temperature is found [49], the parameterizations of the temperature dependence provided in Ref. [49] separately before and after annealing have been used to correct the measurements to the value expected at 0∘*C (in the middle of the Run 2 temperature range). Test-beam measurements taken with ATLAS pixel sensors [3, 51] are also reported in the table but they are more indirect and less precise than TCT measurements.

Both TCT references find that annealing results in a decrease of $\beta$ for electrons and an increase for holes, with a plateau reached after one or two days at 60*∘*C. The two references agree on the values of $\beta$ for electrons, but Ref. [50] finds a smaller $\beta$ value for holes than Ref. [49]. Irradiation with neutrons is found to result in smaller $\beta$-values than irradiation by charged particles [49]. The uncertainties reported in Table 5 do not include an uncertainty of about 10% in the fluence received by the devices, and refer to the average found by fitting measurements made on several devices.

For the simulation results reported in this paper, a value of $\beta_{e}=(4.5\pm 1.5)\times 10^{-16}$ cm2/ns was used for electron trapping and $\beta_{h}=(6.5\pm 1.5)\times 10^{-16}$ cm2/ns for hole trapping. These values were chosen after considering the irradiation conditions of the IBL pixel modules during the LHC Run 2. The range of selected $\beta$ values covers the variation of measurements presented in Table 5 including annealing.

4.6 Ramo potential and induced charge

Charge carrier movement induces a signal on the detector electrodes. The instantaneous current $i$ induced in an electrode by a carrier of charge $q(=\pm e)$ moving at velocity $\vec{v}$ can be calculated by means of the Shockley–Ramo theorem [40, 41]:

[TABLE]

where $\vec{E}_{\text{w}}$ is the Ramo (or weighting) field that describes the induction coupling of the moving charge to a specific electrode. The Ramo field can be calculated by applying a unit potential to the electrode under consideration and zero potential to all other electrodes. Integrating Eq. (12) over a certain drift time, the charge $Q$ induced on the electrode can be expressed as:

[TABLE]

where $\phi_{\text{w}}$ is the Ramo (or weighting) potential with $\vec{E}_{\text{w}}=-\nabla\phi_{\text{w}}$ and $\vec{x}_{\text{start,end}}$ are the initial (final) positions of the charge carrier under consideration. The Ramo potential depends only on geometry and therefore can be computed once prior to any event simulation.

The Ramo potential is calculated using TCAD to solve the Poisson equation; for planar sensors, most of the variation in the Ramo potential is in the $z$ direction, but the $x$ and $y$ dependence must also be included in order to account for charge induced on the neighbouring pixels. Figure 14 shows a slice of the three-dimensional Ramo potential in the centre of the pixel electrode ($y=0$). The vertical line indicates the edge of the pixels: the Ramo potential has sizeable contributions in the neighbouring pixels.

The combination of the Ramo potential and charge trapping is illustrated in Figure 15 for planar sensors. On the electrode of the same pixel in which the electrons and holes originate, the induced charge equals the electron charge if the time to be trapped exceeds the time to drift toward the electrode (Figure 15(a)). The average collected charge is an asymmetric function of the depth inside the sensor because the drift and trapping times are different for electrons and holes and the Ramo potential is very asymmetric: the average fraction is lower far away from the collecting electrode. The charge induced on neighbouring pixels is shown in Figures 15(b) and 15(c). As the trapping time exceeds the drift time, the integrated induced current amounts to the full electron charge in the primary pixel and the charge in the neighbours is zero. For some combinations of starting location and time to trap, the induced charge can have the opposite sign. This happens when holes are trapped very close to the pixel implants.

4.7 3D sensor simulations

In contrast to planar sensors where electrons and holes drift along the depth toward the pixel implant or the back plane, respectively, in 3D sensors charges drift laterally. Columns are etched through the $p$-type silicon bulk (with an initial effective doping concentration of $6\times 10^{11}/\text{cm}^{3}$), and are subsequently either $p^{+}$ or $n^{+}$ doped.161616In the ATLAS IBL, some sensors have columns that extend through the entire bulk while in others they only partially pass through. This section only considers the fully passing through case. In the 3D sensors, two $n^{+}$ columns are shorted so that one pixel corresponds to two $n^{+}$ electrodes for collecting electrons. The $p^{+}$ electrodes are connected to the bias voltage. As the distance between electrodes can be much smaller than the sensor depth, 3D sensors are designed to be more radiation hard than planar sensors due to the reduced drift length (see e.g. Ref. [52] and references therein). This section describes how the digitizer presented in the previous sections can be modified to also accommodate the 3D sensor geometry.

Radiation damage effects for the 3D sensor are implemented in the Perugia model [53] (the Chiochia model cannot be used, as it is designed for $n$-type bulk) with the Synopsys TCAD package [54]. As shown in Figure 16 one-eighth of the sensor is simulated to take advantage of the symmetry within the pixel. In the Perugia model, there are two acceptor traps and one donor trap, with activation energies given by $E_{\text{c}}-0.42$ eV, $E_{\text{c}}-0.46$ eV, and $E_{\text{v}}+0.36$ eV, respectively. The density of traps is predicted to increase linearly with fluence $\Phi$, so for each trap an introduction rate ($g_{\text{int}}$) is defined as: $g_{\text{int}}=N/\Phi$, where $N$ is the trap concentration. In Table 6 the values used for the simulation of 3D sensors reported in this paper are summarized.

Electric field profiles simulated with TCAD are shown in Figure 17. In contrast to planar sensors (Section 4.2.2), the field is nearly independent of $z$ and depends strongly on $x$ and $y$. Therefore, the electric field magnitude is shown as a two-dimensional map for both an unirradiated sensor and a highly irradiated sensor. The $n^{+}$ and $p^{+}$ implants are regions of no field due to their large doping and are modelled as having $0\%$ charge collection efficiency. To illustrate the entire pixel, the one-eighth map that is simulated is tessellated.

The projected time to reach the electrode can be computed analogously to planar sensors (Section 4.3) and is shown in Figure 18. The main difference relative to planar sensors is that electrons and holes follow a non-trivial trajectory due to the more complex electric field. However, a simplification relative to planar sensors is that the electric field is nearly parallel to the magnetic field so the Lorentz angle is negligibly small.

The Ramo potential for 3D sensors is slightly more complex than for planar sensors. In particular, the two $n^{+}$ columns in one pixel are electrically connected and so in the calculation of the Ramo potential, both are held at unit potential while all other electrodes are grounded. Therefore, the calculation requires a relatively large simulation area. This is illustrated in Figure 19. As with the planar sensors, only the immediate neighbours are included in the calculation. The numbers overlaid on Figure 19 show that this is a good approximation.

The various digitizer components are combined in Figure 20 to show the charge induced by electrons as they drift toward the $n^{+}$ electrode. The stochastic path from the same initial electron position near the upper right part of the plot is re-simulated many times. Markers indicate the final location of the electrons (they all have the same initial position). The electrons that travel closer to the electrode before being trapped induce a larger charge (darker markers) than those that are trapped right away.

Since the 3D sensors in the IBL are outside of the acceptance to reconstruct tracks ($|\eta|<2.5$), they have not been studied as thoroughly as their planar counterparts. However, future studies with these sensors will provide an important opportunity to test the 3D digitization model presented here, which has already been used to make projections for the inner tracking detector (ITk) at the HL-LHC [55].

5 Model predictions and validation

5.1 Data and simulation

The models presented in the previous sections are validated by comparing the simulations with data, both in terms of physics predictions. This section presents two key observables for studying radiation damage: the charge collection efficiency and the Lorentz angle. These two quantities are measured as a function of time in Run 2 for the IBL planar sensors. The IBL is well-suited for this test because at the start of Run 2, it was unirradiated. The data were collected in the fall of 2015 and throughout 2016 and 2017. Charged-particle tracks are reconstructed from hits in the pixel detector, silicon strip detector, and transition radiation tracker. Clusters on the innermost pixel layer associated with tracks are considered for further analysis. The IBL is operated with a bias voltage of 80 V, 150 V, or 350 V at a temperature ranging from $-15^{\circ}$C to $+20^{\circ}$C. The analogue threshold is 2550 $e$ with 4 bits of ToT for the digital charge read-out [5]. The ToT is calibrated so that a ToT of 8 corresponds to 16 k$e^{-}$ and a digital threshold of ToT $>1$ is applied. The analogue-to-digital conversion is modelled using the same charge-to-ToT conversion. In practice, the charge-to-ToT conversion is non-linear, especially for low charge. This is most prominent at low charge and therefore required care in interpreting data-to-simulation comparisons in this regime.

Simulated datasets are based on Geant4 [31] with digitization implemented in Allpix [56], which is a lightweight wrapper of Geant4 that is optimized for test-beam analysis and is a powerful test-bench for digitizer development. The figures in this section refer to this set-up as the ‘Stand-alone Simulation.’

5.2 Charge collection efficiency

The collected charge is represented by the most probable value of the charge distribution, which is approximately Landau-distributed [57]. The charge collection efficiency (CCE) is defined here to be the collected charge at one fluence divided by the charge for unirradiated sensors in over-depletion. Figure 21 shows the measured and predicted charge collection efficiencies as a function of integrated luminosity in Run 2. Data points are corrected in order to account for the drift in the ToT calibration. As radiation effects cause the measured ToT to drift with integrated luminosity, regular re-tunings were performed to bring the mean ToT back to the tuning point. This is a feature of the electronics and is not due to a physical change in the charge collection. The drift is approximated as linear and the correction is evaluated in the middle of the run (one period of stable beam) considered. An uncertainty of 30% is assigned to this correction, which is then propagated to the charge collection efficiency value. The size of the correction varies from run to run, but is in general below 5% with a final uncertainty in the CCE of $2$–$3$%. The total uncertainty in the predicted value of the charge collection efficiency is evaluated by taking the squared sum of the differences between the nominal value and the one obtained with the variation of the radiation damage parameters, as explained in detail in Section 4.2.3. An additional uncertainty is due to the trapping constant, as explained in Section 4.5. An uncertainty of 3% [20] is also assigned to the luminosity value of the data points (the horizontal error bars in Figure 21). For the Allpix (‘stand-alone’) simulation points, the integrated luminosity is converted to a fluence using the information presented in Figure 1, and an uncertainty of 15% (see Section 3.1) is assigned to the conversion. As expected, the efficiency drops with luminosity ($\propto$ fluence) due to charge trapping and under-depletion. Part of this loss was recovered after switching the bias voltage in the IBL first from 80 V to 150 V and then later to 350 V, and further increases will be necessary to recover future losses in the charge collection efficiency. A breakdown of the impact of the variations performed to assess systematic uncertainties is reported in Table 7.

5.3 Lorentz angle

The Lorentz angle is determined by performing a fit to the transverse cluster size $F$ as function of the incidence angle of the associated track using the following functional form:

[TABLE]

where $\alpha$ is the incidence angle,171717The angle the tangent vector of the track makes with the vector normal to the sensor surface. $\theta_{\text{L}}$ is the fitted Lorentz angle, $G$ is a Gaussian probability distribution evaluated at $\alpha$ with mean [math] and standard deviation $\sigma$, and $a$ and $b$ are two additional fit parameters related to the depletion depth and the minimum cluster size, respectively. An example input to the fit is shown in Figure 22(a). In general, the simulation does not match the data at very low and high incidence angles, since the simulated points depend on many features of the simulation, but the position of the minimum should depend only on the Lorentz angle. For example, the geometry used for this simulation is simplified and the extreme incidence angles are likely more impacted in the actual geometry. The simulation in Figure 22(a) matches the low incidence angles well, but this is not seen for all fluences; it could be due in part to the uncertainty in the fluence.

The fitted Lorentz angle as a function of integrated luminosity is shown in Figure 22(b). Due to the degradation in the electric field, the mobility and thus the Lorentz angle increase with fluence. This is not true for the Petasecca model, which does not predict regions of low electric field. Charge trapping does not play a significant role in the Lorentz angle prediction. The overall normalisation of the simulation prediction is highly sensitive to the radiation damage model parameters, but the increasing trend is robust. An overall offset (not shown) is consistent with previous studies and appears even without radiation damage (zero fluence) [58], which is why only the difference in the angle is presented.

6 Conclusions and future outlook

This paper presents a digitization model for ATLAS planar and 3D sensors that includes radiation damage effects. Predictions for the fluence at a given integrated luminosity are validated using leakage current data and stand-alone Hamburg model-based calculations. TCAD simulations with effective traps in the silicon bulk are used to model distortions in the electric field caused by exposure to radiation. The impact of annealing is studied by using predictions of the effective doping concentration to adjust the concentration of defect levels in the TCAD simulation. Systematic uncertainties in all aspects of the luminosity-to-fluence conversion and radiation damage model are estimated.

Comparisons between simulations using the radiation damage model and collision data indicate that within the current precision, the fluence-dependence is well-reproduced. The charge collection efficiency gradually decreases with integrated luminosity until there are significant regions of the sensor with a small electric field, at which point there are significant losses. After switching the bias voltage in the IBL from 80 V to 150 V, and then from 150 V to 350 V, the charge collection efficiency significantly increased. The bias voltage will need to be increased further in order to recover the losses in the charge collection efficiency. Both the simulations and the data indicate that the Lorentz angle increases with fluence. The prediction for the Lorentz angle is quite sensitive to variations in the radiation damage parameters, but an increasing trend is a robust prediction of the simulation. With more fluence and annealing, collision data may even be used to further constrain the radiation damage models.

However, the current radiation damage models have known limitations. As an example, in addition to predicting the wrong fluence for space-charge sign inversion, an alternative model (Petasecca) predicts a linear electric field profile that does not qualitatively describe the Lorentz angle dependence on fluence. At the same time, the Hamburg model used to describe annealing does not incorporate effects from a non-uniform space-charge density distribution. Thus far, the Hamburg model provides an excellent description of leakage current data, but this may not hold in the extreme irradiation regime when the electric field profile is sufficiently far from linear. After the next long shutdown currently scheduled to end in 2021, annealing may produce a large enough effect so that the approximate integration of annealing into TCAD presented earlier may no longer be accurate. Another related challenge is that TCAD models such as the Chiochia one need to be improved in order to accommodate an accurate temperature dependence. While future collision data may be used to tune the radiation damage models, additional work may be required to combine the best of the Hamburg and TCAD models. The flexibility of the digitizer model will allow collision data to be used to validate and test new ideas in the future. A version of the digitizer model is publicly available on GitHub [59] and the model is also implemented in the ATLAS software framework (ATHENA) in releases designed for simulating the data in 2016 and beyond. The ATLAS code can be found on GitLab [60].

Even though the innermost pixel layers have only been exposed to $\Phi\approx\times 10^{15}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$, radiation damage effects are already measurable. The projected fluence on the IBL at the end of LHC operation (300 fb*-1*) is about $\Phi=2\times 10^{15}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$; the sensors on the innermost layer of the upgraded ATLAS tracker (ITk) will have to withstand $\Phi\sim 1\times 10^{16}$ $\text{n}_{\text{eq}}/\text{cm}^{2}$ [61]. For such high fluences, pixel sensor radiation damage will be an important aspect of operations (setting the high voltage, deciding the time spent warm, etc.) and track reconstruction. The simulation framework presented here can be used to inform both online and offline performance for the current pixel detector as well as for making important design decisions for the upgraded ATLAS detector that must survive the harsh HL-LHC radiation environment.

Acknowledgements

We thank CERN for the very successful operation of the LHC, as well as the support staff from our institutions without whom ATLAS could not be operated efficiently.

We acknowledge the support of ANPCyT, Argentina; YerPhI, Armenia; ARC, Australia; BMWFW and FWF, Austria; ANAS, Azerbaijan; SSTC, Belarus; CNPq and FAPESP, Brazil; NSERC, NRC and CFI, Canada; CERN; CONICYT, Chile; CAS, MOST and NSFC, China; COLCIENCIAS, Colombia; MSMT CR, MPO CR and VSC CR, Czech Republic; DNRF and DNSRC, Denmark; IN2P3-CNRS, CEA-DRF/IRFU, France; SRNSFG, Georgia; BMBF, HGF, and MPG, Germany; GSRT, Greece; RGC, Hong Kong SAR, China; ISF and Benoziyo Center, Israel; INFN, Italy; MEXT and JSPS, Japan; CNRST, Morocco; NWO, Netherlands; RCN, Norway; MNiSW and NCN, Poland; FCT, Portugal; MNE/IFA, Romania; MES of Russia and NRC KI, Russian Federation; JINR; MESTD, Serbia; MSSR, Slovakia; ARRS and MIZŠ, Slovenia; DST/NRF, South Africa; MINECO, Spain; SRC and Wallenberg Foundation, Sweden; SERI, SNSF and Cantons of Bern and Geneva, Switzerland; MOST, Taiwan; TAEK, Turkey; STFC, United Kingdom; DOE and NSF, United States of America. In addition, individual groups and members have received support from BCKDF, CANARIE, CRC and Compute Canada, Canada; COST, ERC, ERDF, Horizon 2020, and Marie Skłodowska-Curie Actions, European Union; Investissements d’ Avenir Labex and Idex, ANR, France; DFG and AvH Foundation, Germany; Herakleitos, Thales and Aristeia programmes co-financed by EU-ESF and the Greek NSRF, Greece; BSF-NSF and GIF, Israel; CERCA Programme Generalitat de Catalunya, Spain; The Royal Society and Leverhulme Trust, United Kingdom.

The crucial computing support from all WLCG partners is acknowledged gratefully, in particular from CERN, the ATLAS Tier-1 facilities at TRIUMF (Canada), NDGF (Denmark, Norway, Sweden), CC-IN2P3 (France), KIT/GridKA (Germany), INFN-CNAF (Italy), NL-T1 (Netherlands), PIC (Spain), ASGC (Taiwan), RAL (UK) and BNL (USA), the Tier-2 facilities worldwide and large non-WLCG resource providers. Major contributors of computing resources are listed in Ref. [62].