Microscopic origin of the acceptor removal in neutron-irradiated Si detectors - An atomistic simulation study

.


Introduction
Si sensors fabricated on lowly doped p -type substrates are a technological solution adopted for the next generation of particle detectors. P -type sensors show an improved radiation hardness that makes them suitable to withstand the harsh radiation fields expected in high energy physics experiments [1] . One of the main effects of irradiation in Si devices is the so-called displacement damage. It arises from the generation of Si self-interstitials (I's) and vacancies (V's) when Si atoms are moved from their lattice positions. Radiation-induced defects have a strong impact on the detector electrical behavior due to an effective dopant concentration (N eff ) whose value may notably differ from that before irradiation [2] . As a consequence, the resistivity of the detector depleted region changes and the net space charge is modified, affecting the electric field profile and the full depletion voltage (V FD ). Radiation-induced changes on N eff represent a serious challenge * Corresponding author.
Recent works have addressed the experimental characterization of N eff upon hadron irradiation in p -type sensors by measuring the variations on resistivity and V FD with fluence [ 3 , 5 , 7-9 ]. As many detectors are expected to operate in cooled conditions for a lifetime of several years, accelerated annealing tests are used to analyze system dynamics. It is commonly assumed that a stable space charge, i.e. a constant value of N eff over time, can be obtained after annealing for 80 min at 60 °C [10] . In these conditions, N eff is strongly dependent on irradiation fluence. At low fluences, N eff follows an exponential decay with fluence attributed to an acceptor removal process [ 2 , 10 ], i.e. the deactivation of dopants (usually B atoms) by displacement damage. The fast decay of N eff is characterized by the so-called "acceptor removal constant" that reduces as the initial dopant concentration (N eff,0 ) increases [4] . The mechanisms behind this dependence are still unknown, and it remains as one of the most intriguing features of acceptor removal. At high fluences, N eff increases linearly with irradiation fluence, which is associated to the formation of defect-related deep acceptor levels [ 10 , 11 ]. Electrically active defects can be experimentally characterized to some extent by Deep Level Transient Spectroscopy (DLTS) and Thermally Stimulated Current (TSC) techniques [12] . However, the large variety of defect sizes and topologies hinders the identification of the defects responsible for the spectra peaks, and their usually small size prevents the direct observation by microscopy.
The modeling of radiation-induced changes on N eff is complicated because of the large amount of defect-related energy levels observed after irradiation, and the difficulty in linking macroscopic effects with microscopic defects [ 12 , 13 ]. Most models attempt to simulate the electrical behavior of irradiated devices by posing a phenomenological description of the process. The acceptor removal is usually modeled with a few first order equations for defect reactions with dopants and impurities, and the use of several parameters to account for the generation of displacement damage and the formation of B defects [ 3 , 14 , 15 ]. The increase of N eff at moderate fluences is simplified by assuming that the effects of displacement damage can be effectively represented by a few donor/acceptor energy levels in the Si bandgap, with parameters fitted to reproduce the macroscopic degradation of the detector electrical properties [ 16 , 17 ]. Although available empirical models could be fitted to reproduce the macroscopic behavior of devices to some extent, they do not provide insight on the physical mechanisms underlying the changes of N eff with fluence. Detailed multi-scale simulation schemes have been proposed to analyze the effects of irradiation induced defects on the device electrical and optical characteristics. Nevertheless, previous works focused on residual intrinsic defects and did not include dopants [18][19][20] . We presented a preceding atomistic model to analyze B deactivation in irradiated Si [21] , but a comprehensive understanding of the acceptor removal process is still lacking.
In this work, we address the acceptor removal process by atomistic simulations of the generation of displacement damage under neutron irradiation, and the interactions between defects and dopants. We analyze the microscopic origin of the acceptor removal mechanism, provide an explanation for the dependence of the acceptor removal constant on initial resistivity, and give some clues to identify those intrinsic defects that could act as deep acceptors.

Simulation scheme
We perform atomistic simulations to analyze the displacement damage resulting from neutron irradiation on p -type Si substrates, and its effect on the dopant profile. The energy spectrum of reactor neutrons typically spans from 10 −8 to 10 MeV [22] but, in our simulations, we consider monoenergetic neutrons of 1 MeV as it is the standard reference for the Non-Ionizing Energy Loss (NIEL) scaling, and it is close to the NIEL weighted mean energy of the Jožef Stefan Institute (JSI) reactor (1.7 MeV) [13] . A flux of 10 12 n cm −2 s −1 is assumed, which is in the range of those commonly used at the JSI reactor [22] . P -type substrate doping is simulated by considering B atoms randomly distributed, with five different concentrations ranging from 10 13 to 10 17 cm −3 . Cubic simulation cells of at least 3 × 3 × 3 μm 3 are used to get enough statistics for all B concentrations. Si recoils are generated at random within the simulation volume, and periodic boundary conditions are applied in all directions to analyze the kinetics in the bulk. Irradiation is performed at RT (25 °C), covering 1-MeV neutron equivalent (n eq ) fluences from 10 12 to 10 16 cm −2 for a time ranging from 1 to 10 4 s, respectively. After each irradiation step the sample is annealed at 60 °C for 80 min, following the so-called CERN scenario measurement technique [23] .
Several tools are combined in a multi-scale scheme so that the effects of neutron irradiation on p -type Si can be analyzed at an atomic level and, at the same time, macroscopic dimensions and time scales can be reached. The energy distribution of primary knock-on atoms (PKA) due to the neutron irradiation is calculated with the SPECTRA-PKA code [24] using the TENDL-2017 nuclear data library. Monoenergetic 1 MeV neutrons result in a total recoil matrix with a PKA energy spectrum that spans up to 150 keV with an average energy of 48.5 keV. Simulated neutron fluences correspond to PKA fluences from 6 × 10 7 to 6 × 10 11 cm −2 , which implies that in our simulation cell and for the highest fluence more than 54,0 0 0 PKAs are simulated. Following the PKA energy spectrum, the corresponding number of recoils for each PKA energy is simulated with MARLOWE, a tool based on the binary collision approximation that provides the coordinates of the Si I's and V's [25] . The positions of Si I's and V's are transferred to the object kinetic Monte Carlo (kMC) process simulator DADOS, which only simulates the kinetics of defects and dopants during annealing and ignores the lattice [26] . In the kMC technique, system dynamics is modeled by performing random events whose probability is determined by the energetics of the defects and dopants involved in the process.
The insight gained into the mechanisms of B deactivation by implantation damage and their implementation in kMC models as the result of an extensive research on junction formation in Si devices [27] , can be applied to the system analyzed in this work. However, while the fabrication of B-doped junctions often involves high local concentrations of defects and dopants and intense thermal budgets, in this study damage is diluted, B concentration is low and the annealing is reduced. Under the low thermal budget used in the aging annealing of detectors, system dynamics is governed by low activation energy events, mainly related to the diffusion of mobile species. The main interactions between defects, dopants and impurities considered in the simulations of this work are summarized in Table 1 , and their energetics are provided in Tables S1-S5 in Supplementary Material. The interaction of mobile Si I's and V's is assumed to cause their mutual annihilation. When defects of the same type interact they form agglomerates of Si I's or V's (called clusters) that are mostly considered as immobile (only the di-interstitial (I 2 ) and tri-interstitial (I 3 ) defects can diffuse).
A boron-interstitial (B i ) is formed when a Si interstitial interacts with a substitutional B atom. B i can also diffuse and interact with other defects and dopants resulting in the formation of boroninterstitial clusters. They are referred as B n I m , where n and m denote the number of Si I's and B atoms, respectively, and can grow in size by incorporating additional Si I's or B atoms [28] . We assume that substitutional B atoms are electrically active, contributing as shallow acceptors, while B i and B n I m are considered electrically inactive and cause B deactivation. In a previous work, we explored the capabilities of our multi-scale approach to describe the damage induced by neutron irradiation, and analyzed the acceptor removal for a unique B concentration [21] . We showed that the diffusion of the I 2 and I 3 defects, which is ignored in B junction formation, plays a role on B deactivation at room or slightly elevated temperatures.
C and O impurities are incorporated in a basic model through first-order reactions. They act as traps for Si I's, V's and B i leading to the formation of C i , VO and B i O defects, but higher order reactions are not included. We assume that B i O complexes also deactivate B atoms, but they are not assigned any explicit carrier contribution. Therefore, in our code B i , B i O and B n I m are considered electrically inactive defects and they all contribute to B deactivation. Simulations were done considering 10 17 O cm −3 and 2 × 10 15 C cm −3 , in the range of impurity concentration in standard Si wafers [ 11 , 13 ], for all B concentrations and neutron fluences indicated above. Although more detailed models for the interactions of C and O with intrinsic defects have been proposed, for the conditions analyzed in this work our simplified model correctly describes the amount of Si I's and V's immobilized by C and Interstitial-rich C complexes, such as CI 2 or CI 3 , are shown to be energetically unfavorable [29] , and C-rich complexes are unlikely to form for the low C concentration simulated. Ab initio calculations suggest that O can bind to vacancy clusters leading to stable V n O m complexes [30] . However, this will not alter the overall number of immobilized vacancies, as V's clusters are already stable for the temperature range of our simulations. As a benchmark for our model, we compare the experimental introduction rates for VO derived from DLTS with those in our simulations. Depending on experimental conditions, introduction rates of 0.1, 0.63 and 3.4 cm −1 have been reported [31][32][33] , usually for low fluence irradiation. In our simulations with O and C (10 15 B cm −3 , 10 17 O cm −3 , 2 ×10 15 C cm −3 ) irradiated with 10 12 n eq cm −2 we obtain a VO concentration of 1.3 × 10 12 cm −3 , which corresponds to an introduction rate of ∼1.3 cm −1 . This value is consistent with experimental results, although the amount of VO in simulations may be overestimated, as some of them may evolve to higher order V n O m complexes that are not considered in our model.
Unlike empirical models that use first order equations for defect-dopant reactions and fitting parameters for the amount of displacement damage and the probability of defect-dopant interactions, in our simulation scheme, the amount, distribution and topology of damage and defect-dopant interactions are naturally derived from atomistic simulations, with no previous assumptions. Our model takes into account the complex interactions between Si I's and B atoms, including higher order reactions, which allows a detailed analysis of B deactivation processes. However, our simulation scheme has some limitations. The kMC simulator does not explicitly include the charge state of dopants and defects. For the discussion on the electrical activity of defects we assume those states obtained by ab initio calculations reported in literature. Impurities (C, O) are described with a simplified model that captures their role as traps but not the complex scenario that may arise from the interactions between intrinsic defects, dopants and impurities (C n I m , V n O m , C-O, B-O, B-C complexes).

Results and discussion
We have analyzed the acceptor removal process by simulating B deactivation under neutron irradiation. We use experimental data available in literature [ 7 , 8 , 34 ] and our own experimental results for a radiation-hard power switch described elsewhere [ 21 , 35 ] to validate our model [21] . Estimated N eff values from the measured drain currents in our switch, along with other experimental N eff obtained from V FD measurements of irradiated devices [ 7 , 8 , 34 ] are plotted in Fig. 1 . At first, N eff significantly reduces with fluence, but at high fluences a linear increase is observed. This increase has been associated with the contribution to space charge of electrically active defects, which have a net acceptor behavior, increasing V FD . The evolution of the experimental N eff in our switch is slightly different. It presents a small plateau at low fluences, then a steeper decay with fluence and a final saturation at high fluences, but no increase is observed. In our case, N eff is not derived from V FD measurements that detect the amount of space charge, but from bulk resistivity values, which is actually a measure of the concentration of free carriers. At low fluences, the concentration of free carriers equals the space charge concentration (since space charge is mostly due to ionized dopants), but at high fluences the contribution of electrical active defects to space charge exceeds that of dopants. The simulated N eff (open diamonds in Fig. 1 ) corresponds to the concentration of substitutional B atoms, i.e. B atoms that are still active and can contribute as acceptors. At low fluences, our simulated data follow the same trend as experimental N eff values. To get a global picture of the evolution of N eff with fluence, we have added to the simulated N eff values the contribution of deep acceptors (dashed line) (calculated as the experimental introduction rate 0.02 cm −1 × fluence [ 10 , 11 ]). The resulting N eff curve (plotted as solid diamonds) provides a consistent description of the process in the whole fluence range.
The dependence of N eff on fluence is usually parametrized following the so-called "Hamburg model" [36] as where φ eq is the 1-MeV NIEL equivalent fluence, N eff,0 the initial acceptor concentration, N C and N C,0 the concentration of removed acceptors and removable acceptors (N C,0 ≤ N eff,0 ), respectively, c A the acceptor removal constant, and g c the introduction rate of stable defect-related deep acceptors. The Hamburg model takes into account both the sharp reduction of N eff with fluence at low fluences, and its linear raise at higher fluences. The acceptor removal process can be better analyzed by defining the parameter φ A as the inverse of the removal constant ( φ A = 1/c A ). φ A represents the fluence at which N C has reached 63% of its final value (N C,0 ) and identifies the steep-decay region in N eff curves. It can be envisioned as a figure of merit of the detector radiation hardness as it refers to the fluence that the detector can endure before significant dopant deactivation occurs. In the following sections, we present simu-   [ 6-9 , 38 , 39 ], LGAD [ 3 , 37 ], sensors made from EPI Si [ 34 , 37 ] and Pad diodes [39] .
lations of acceptor removal under neutron irradiation for a wide range of substrate resistivities. We use experimental data from literature to analyze the physical mechanisms that explain the values of φ A and g c , and the dependence of φ A on N eff,0 , as these are the key parameters that determine the evolution of N eff with fluence.

Acceptor removal parameter φ A
The acceptor removal constant c A has been systematically reported to decrease with N eff,0 (i.e. φ A increases with N eff,0 ) [ 3-5 , 9 , 37 ]. To analyze this dependence, we have performed a total number of 90 simulations, including 9 neutron fluences and 5 N eff,0 values, with and without impurities (O, C). Dopant concentration ranges between 10 13 and 10 17 B cm −3 , covering the high-and lowresistivity substrates of standard Si sensors and novel detectors (LGAD, HV-CMOS), respectively. In Fig. 2 we show the simulated relative effective dopant concentration, i.e. the ratio N eff /N eff,0 , as a function of 1-MeV n eq fluence for the doping range mentioned before. Our results reveal that as N eff,0 increases the sharp decay region is shifted to larger fluences, consistently with experimental observations. For instance, for 10 13 B cm −3 N eff is reduced to a 50% of its original value after a fluence of 5 × 10 13 n eq cm −2 , while for 10 17 B cm −3 a fluence of 10 16 n eq cm −2 is required to achieve the same reduction in N eff /N eff,0 . Simulated N eff curves are fitted to Eq. (1) using N C,0 and φ A as fitting parameters (solid lines in Fig. 2 ). The extracted values of φ A as a function of N eff,0 are compared in Fig. 3 to those obtained from the experimental removal constants c A ( φ A = 1/c A ) available in literature. Experimental data correspond to the analysis of N eff obtained from V FD measure- ments, and include different detector architectures as HV-CMOS sensors [ 6-9 , 38 , 39 ], LGAD [ 3 , 37 ], sensors made from EPI Si [ 34 , 37 ] and Pad diodes [39] . φ A values obtained from simulations are comparable to the experimental ones and follow a sublinear trend with N eff,0 according to equation φ A = 2 × 10 6 N eff,0 0.56 cm −2 .

Role of O and C on acceptor removal
The presence of O and C in Si wafers in concentrations comparable or higher than N eff,0 [ 11 , 13 ] may play a role on acceptor removal since they interact with defects and dopants. It is well established that the existence of O leads to the formation of VO and higher order V-O complexes by capturing V's [40] , while C atoms react mainly with Si I's via the formation of C i and C-I complexes [41] . Several studies indicate that B i has a strong affinity for pairing with O, resulting in B i O complexes [ 42 , 43 ]. Traditionally, B i O is assumed to behave as a donor in the upper part of the Si bandgap and is assigned a key role on acceptor removal [ 42 , 43 ]. However, recent studies attribute tentatively this donor level to B-I defects [44] , and the large scattering of B i O electrical characterization suggests a more complex scenario in which higher order defects involving B, O and Si I's may exist [45] . To analyze the role of impurities on the removal process new simulations were done assuming 10 17 O cm −3 and 2 ×10 15 C cm −3 , consistent with impurity concentrations in standard Si substrates. The resulting N eff /N eff,0 and φ A values are also depicted in Fig. 2 and Fig. 3 (open symbols), respectively. The introduction of these concentrations of O and C has a small effect on the amount of removed acceptors, and simulations with and without impurities lead to similar removal parameters.
In order to elucidate B deactivation mechanisms, we have quantified in Fig. 4 the fraction of B ([B]/N eff,0 ) that remains substitutional (i.e. active B), in samples with and without impurities after 10 15 n eq cm −2 . We also show the relative weight on the total amount of deactivated B ([B]/N C ) of B i , B n I m and B i O complexes. For low N eff,0 , in the sample without impurities, deactivated B is almost equally distributed between B i and B n I m . However, in the sample with O and C, most deactivated B belongs to B i O complexes, the amount of B n I m has notably reduced and there is almost no B i . At high N eff,0 , the large amount of B atoms enhances B n I m formation becoming the main mechanism of acceptor removal in both samples.
A detailed analysis of our simulations indicates that in the sample without O and C, B deactivation starts with the formation of B i and small B n I m (BI 2 , BI 3 ) through the interaction of mobile Si I's defects (I, I 2 , I 3 ) with B (see Table 1 ). As N eff,0 increases, diffusing B i defects are likely to interact with other B atoms (substitutional B, B i , B n I m ). Thus, the fraction of B in B n I m increases at the expense of B i . In the presence of O and C, for low N eff,0 , oxygen concentration is significantly higher than that of B, which enhances the interaction of B i with O forming immobile B i O complexes. This notably reduces B i population and limits the growth path for larger B n I m , as in our model the interaction of Si I's or B i with B i O complexes is not considered. Besides, a high concentration of VO defects is also formed that can annihilate mobile Si I's (I mob ), hindering the formation of small B n I m through the diffusion of I 2 and I 3 . Regarding C atoms, their main role is the trapping of Si I's in C i reducing B diffusion and clustering [3] . However, in our sample the effect of C is mitigated by the relatively low C concentration and a much higher O concentration (usual in standard Si wafers). As many V's are immobilized at VO defects, I-V annihilation is reduced. This increases the amount of available mobile Si I's and compensates the trapping of Si I's in C i . Therefore

Dependence of removal parameter on substrate resistivity
Experimental and simulated results reported in Fig. 3 show that as N eff,0 increases larger neutron fluences (higher values of φ A ) are required to deactivate the dopants. Previous attempts to explain the acceptor removal process based on defect kinetics did not provide a fully satisfactory description of experimental data [ 14 , 15 ], and this dependence stills remains a riddle. In order to shed light on this issue we have studied the atomistic mechanisms behind dopant deactivation.
From our simulations, we analyze the damage resulting from 1 MeV neutron irradiation and the effect of N eff,0 and fluence on the deactivation process. In Fig. 5 we show the concentration of displacement damage (Si I's and V's) generated in irradiation cascades (dashed line) as provided by Marlowe before any annealing is done. The concentration of I mob (circles and solid line) and Si I's and V's clusters (squares and dotted line) as a function of fluence is also plotted. I mob and cluster concentrations are computed after irradiation steps and are averaged values from simulations with 10 13 to 10 17 B cm −3 , although standard deviation (error bars in the figure) is negligible. The generated displacement damage increases linearly with fluence as it is only determined by the num-ber and energy of PKAs. The great difference between the number of generated defects and those remaining (I mob or clusters) is a clear indication that during RT irradiation intense I-V annihilation has occurred. The amount of I mob and clusters is similar and both increase sublinearly with fluence ([I mob ] = 8227 φ eq 0.77 cm −3 ; [Clusters] = 145 φ eq 0.88 cm −3 ) due to an enhanced annihilation at high fluences as generated defects are in close proximity. The concentration of removed acceptors N C (diamonds) and the acceptor removal parameter φ A obtained for each N eff,0 (open triangles) are also plotted in the figure. We can see that N C initially increases with fluence but it saturates at high fluences when most dopants have already been deactivated (as reported in Fig. 2 ). The saturation onset is shifted to higher fluences as N eff,0 increases and is consistent with the extracted φ A values.
The acceptor removal is caused by the interaction of I mob with substitutional B atoms but V's and Si I's clusters formed along the cascade can annihilate or trap I mob , respectively, hindering them from reaching B atoms. Consequently, the efficiency of the removal process relies on two parameters: the concentration of mobile Si I's remaining after irradiation ([I mob ] ), and the probability of these Si At low fluences and high N eff,0 the amount of B atoms notably exceeds that of annihilating or trapping defects ([B] >> [V cl + I cl ]) and ρ I-B approaches unity. Consequently, the number of I mob,eff , and therefore the amount of B that can be potentially deactivated, is determined by the number of available mobile Si I's (N C ∼ I mob,eff ∼ I mob ). This corresponds in our simulations ( Fig. 5 ) to the region where N C curves (for different N eff,0 ) grow with fluence with a similar slope as I mob , and the values of N C and I mob are comparable. In this regime, the number of I mob is low due to the intense intracascade I-V annihilation, and dopant deactivation is limited by the availability of diffusing Si I's. On the contrary, at large fluences and low N eff,0 , the amount of defect clusters is much  ∝ [B]). This corresponds to the saturated region in Fig. 5 in which N C is described by horizontal lines whose value is determined by N eff,0 . In this situation, there are enough I mob but the presence of a large amount of defect clusters reduces the probability of their interaction with B. The onset of the saturation regime occurs when the amount of I mob and defect clusters is comparable to that of B ([B] ∼ I mob ∼ [V cl + I cl ]). The acceptor removal parameter φ A can be regarded as the fluence required to produce a number of I mob similar to N eff,0 , resulting in significant dopant deactivation and the onset of saturation. Thus, φ A increases with N eff,0 because more I mob are required to deactivate the dopants.
In Fig. 5 we observe some deviation from the previous simplified description. On one hand, for low fluences and dopant concentrations, N C lies below the I mob line and the curves approach I mob as N eff,0 increases. This is mainly due to the inhomogeneity of neutron irradiation damage: defects are close to each other along the cascade trajectory but there is a large distance between neighboring cascades, especially at low fluences. Within 50 keV recoils, each defect (Si interstitial or vacancy) has on average 4.6 neighboring defects within a radius of 2 nm, which is much shorter than the mean inter-cascade distance (403 nm for 10 13 n cm −2 corresponding to a PKA fluence of 6 × 10 8 cm −2 ). As a result, intracascade defect annihilation and clustering is favored, compared to what it would be obtained if the same V's and Si I's concentration were homogeneously distributed [46] . It reduces the actual probability of I mob reaching B atoms and, therefore, the efficiency of the removal process. This effect is made up by increasing the amount of I mob required to deactivate the dopants, which shifts φ A to higher fluences than expected. On the other hand, at high fluences and dopant concentrations N C values slightly exceed the number of I mob . It indicates the existence of more efficient deactivation mechanisms in which each I mob can deactivate several B atoms (I + B → B i ; B i + B → B 2 I). In this case, the high B concentration favors the formation of B n I m with more than one B atom which boosts the removal process. This mechanism reduces the number of required I mob and shifts φ A to lower fluences. As a consequence, although φ A was expected to follow a similar trend as I mob , the inhomogeneity of damage at low fluences and the formation of B-rich B n I m at high N eff,0 result in the sublinear dependence of φ A on N eff,0 obtained in our atomistic simulations and reported in the experiments gathered in Fig. 3 .
Although the overall effect of the simulated concentrations of O and C (when [O] >> [C]) on acceptor removal is limited, they are responsible for slight modifications on defect dynamics. In our simulated samples enriched with O and C, the formation of VO and C i by trapping some V's and Si I's initially prevents them from annihilating and increases the amount of mobile V's and Si I's after irradiation. However, during annealing VO and C atoms can annihilate or trap I mob , limiting its probability of interacting with B. These two effects compensate to a large extent since there are initially more I mob but they are more likely to annihilate in a later stage. It explains why in our model, and for the impurity concentrations used in this work, the amount of deactived B is hardly affected by the presence of O and C.

Microscopic origin of deep acceptors
Experiments reveal that at large irradiation fluences N eff increases almost linearly with fluence, following an acceptor generation rate of g c ∼ 0.02 cm −1 [ 10 , 11 ], as indicated in Eq. (1) . N eff may even exceed N eff,0 , which suggests that other defects in addition to dopants are contributing with negative space charge. Besides, similar values of g c have been reported under neutron irradiation independently of substrate resistivity, in both p -and n -type substrates, and for materials with different O contents [47] . Consequently, the increase of N eff with fluence has been attributed to the formation of radiation-induced electrically active intrinsic defects, with a joint contribution of several donor and deep acceptor levels, the latter ones being predominant [ 12 , 13 ]. This supports the idea that in neutron-irradiated p -type Si N eff is determined by two competing mechanisms: an initial acceptor removal at low and medium fluences caused by B deactivation by Si I's, and the formation of deep acceptors at large fluences associated to the accumulation of intrinsic defects. Extensive research has been devoted to identify the defects most plausible candidates to behave as deep acceptors. Defects proposed from experimental and theoretical studies are in general vacancy-like, and include small vacancy defects or oxygenrelated defects [ 33 , 48 ], clustered or extended defects [ 12 , 13 ] and small amorphous regions [49] . The linking between microscopic defects and macroscopic effects is difficult because a large variety of irradiation-induced defects with different sizes and configurations simultaneously coexists, and they may even contain dopants and impurities. We combine the information obtained from the experimental characterization, ab initio calculations and our kMC simulations to get insight into the main contributors to N eff as deep acceptors.
Experimental defect characterization relies on the use of different irradiation particles to bias damage towards the predominance of point defects or clusters, the controlled introduction of impurities (mainly C and O), and the annealing of irradiated samples to determine the thermal stability of defects. DLTS and TSC techniques characterize electrically active defects, providing parameters such as the associated energy level in the Si bandgap, the donoror acceptor-like behavior, the capture cross-section for carriers and the estimated concentration. However, the identification of specific defects responsible for the peaks of the DLTS/TSC spectrum is not straightforward. Three main acceptor levels were detected, labeled as H116K, H140K y H152K, with energy levels in the lower part of the Si bandgap and a generation rate of ∼ 0.04 cm −1 [13] . These acceptor levels were associated to bistable defect clusters, presenting two configurations with different electrical activity. As their behavior resembled that of small V's clusters these levels were tentatively attributed to clusters with more than 3 V's [12] .
Ab initio simulations allow the characterization of small defect clusters with a previously defined configuration, providing their formation energy and the energy levels introduced in the Si bandgap. Nevertheless, the underestimation of the Si bandgap results in some uncertainty on the calculated energy levels [50] , and the high computational cost limits the electrical behavior studies to small sizes: up to the tetra-interstitial and the hexa-vacancy. Usually, donor levels are obtained for Si I's clusters [ 51 , 52 ], and acceptors levels in the upper part of the Si bandgap (some close to mid-gap) for V's clusters [50] . To our knowledge, there is no theoretical evidence of specific defects with deep acceptor levels consistent with those of H116K, H140K and H152K peaks. The energy levels of small V's clusters derived from ab initio calculations are compiled in Table 2 [ 50 , 53 , 54 ]. Vacancy clusters larger than V 2 are known to be bistable defects existing in two configurations: Partof-Hexagonal-Ring (PHR) and Four-Fold Coordinated (FFC), FFC being the most stable one [50] . For V 3 both configurations have been observed: PHR V 3 is detected at RT, while FFC V 3 is formed after annealing and presents an acceptor level at 75 meV below conduction band edge ( E c ) [54] . The deepest acceptor levels for V's clusters reported from ab initio calculations are 0.46 eV (PHR V 3 −/0 ), 0.54 eV (FFC V 4 = /and V 4 −/0 ), 0.47 eV (FFC V 5 = / − ) and 0.45 eV (FFC V 5 −/0 ) below E c . Note that an acceptor level at E c -0.7 eV was obtained for PHR V 4 −/0 in Ref. [50] , which lies within the lower part of the bandgap, but the authors point out that it may be overestimated in their calculations.
Our kMC simulations provide the concentration and size distribution of Si I's and V's clusters as the result of defects diffusion and interactions in neutron-irradiated p -type Si. Electrically active defects that can significantly contribute to N eff should be present in a concentration similar or higher than that corresponding to the experimental factor g c . In Fig. 6 we report the concentration versus fluence of the most abundant V's and Si I's clusters in a simulation with N eff,0 = 10 15 B cm −3 , and the corresponding fitted power equations as dashed lines (fitting parameters are compiled in Table S6 and Fig. S1 in Supplementary Material). The solid black line corresponds to the effective deep acceptor concentration estimated from the experimental generation rate g c ∼ 0.02 cm −1 . It represents the net contribution of donor and deep acceptor levels, while a higher introduction rate of ∼ 0.04 cm −1 was obtained for the main acceptor levels reported in DLTS/TSC [13] . The figure shows that small clusters are predominant and as cluster size increases their concentration is notably reduced. To estimate whether the amount of clusters we get agrees with experiments, we compare our simulation results against the V 2 concentration obtained from DLTS measurements. Introduction rates of 4.7 cm −1 for V 2 [33] , and 0.3 and 1.4 cm −1 for V 2 ( = ,-) and V 2 (-,0) , respectively [55] , have Table 2 Energy level (E, below the conduction band edge), occupancy ( η) and associated introduction rate of negative space charge (g cl ) of small vacancy clusters considering σ e = 2 × 10 −15 cm 2 and σ h / σ e ratios of 10 and 100. Energy levels are extracted from ab initio calculations available in literature [ 50 , 53 , 54 ]. η lower than 0.01% and g cl lower than 0.001 cm −1 are denoted by "−".  been reported. The V 2 concentration shown in Fig. 6 after 10 12 n cm −2 is 1.8 × 10 12 cm −3 , which corresponds to an introduction rate of 1.8 cm −1 , in the order of those derived from DLTS at low irradiation fluences. According to our simulations, only small clusters with a size up to 5-6 defects are abundant enough compared to the deep acceptor concentration deduced from the experimental generation rates. Intrinsic defects that play a major role on the increase of N eff with fluence must exist in a concentration high enough, and have acceptor energy levels in the lower part of the bandgap or close to mid-gap to have a significant occupancy. We estimate the contribution of an electrically active defect to N eff from the concentration obtained by kMC simulations and its steady-state occupancy. This is calculated from the Shockley-Read-Hall statistics as done in Ref. [13] , assuming full depletion (zero free carriers). We use energy levels in the bandgap reported by ab initio calculations and electron ( σ e ) and hole ( σ h ) capture cross sections derived from DLTS/TSC. For acceptors in the upper part of the bandgap, σ e in the order of 2-4 × 10 −15 cm 2 has been reported [12] . Continuum models use σ h as a fitting parameter assuming σ h / σ e ≈ 10 [16] as few experimental values are known (only σ h / σ e ≈ 100 for PHR V 3 −/0 ) [12] . In Table 2 we show the energy levels (E, below E c ) reported in literature for small vacancy clusters (V 2 -V 6 ) with different charge states, and the calculated occupancy ( η) assuming σ e = 2 × 10 −15 cm 2 and using two different σ h / σ e ratios of 10 and 100. For each cluster size, the total defect concentration is distributed among the different char ge states including the neutral state (not shown in Table 2 ). FFC V 4 clusters show a high occupancy as their energy levels are close to mid-gap ( E c -0.54 eV). Considering cluster concentrations reported in Fig. 6 and the obtained occupancy, we determine the introduction rate of negative space charge associated to each cluster (g cl ) (singly and doubly charged defects contribute with one or two charges to N eff , respectively). Our results indicate that, in full depletion, the contribution of FFC V 4 = / − is dominant and we even overestimate the increase of N eff with irradiation fluence. Other small V's clusters (FFC V 4 −/0 , PHR V 3 −/0 , and FFC V 5 −/0 ) contribute to a lower extent and some compensation with donor-like defects may also occur. Although we have not found a correspondence with the acceptor levels H116K, H140K, H152K measured by DLTS/TSC, our analysis is consistent with empirical radiation damage models of irradiated p -type sensors, which require acceptor levels in the range of 0.42-0.52 eV below E c (similar to the energy levels of V 3 , V 4 , V 5 clusters) to fit their experimental electrical behavior [37] .

Conclusions
We present the first TCAD study that addresses the acceptor removal process in neutron-irradiated p -type Si by means of atomistic simulations of defect-dopant interactions. The acceptor removal can be successfully modeled by simulating B deactivation induced by Si I's generated by irradiation. We combine several simulation techniques to describe neutron irradiation damage and dopant deactivation processes without specific fitting parameters or previous assumptions. Our defect-dopant atomistic model, widely tested in Si process simulation, is extended to include I 2 and I 3 as mobile defects as they play a key role on B deactivation under the operational conditions of particle detectors.
From the analysis of N eff on a wide range of N eff,0 (10 13 -10 17 B cm −3 ) we obtain acceptor removal parameters that are consistent with those experimentally reported and with the same dependence on substrate resistivity. Our study reveals that the efficiency of the acceptor removal process strongly depends on the relative population of I mob , defect clusters and B atoms. We show that the increase of the acceptor removal parameter φ A with N eff,0 is due to the limited amount of I mob at low fluences and the reduced probability of I-B interaction as fluence increases. The shortage of I mob is caused by intense I-V annihilation within irradiation cascades, and by the trapping of Si I's in clusters that prevents them from interacting with dopants. Defect clusters play a key role by annihilating and trapping I mob since, contrary to impurities and other traps for Si I's whose number remains constant during irradiation, cluster concentration increases with fluence and may even exceeds that of I mob . To provide a detailed description of the acceptor removal mechanism, B deactivation models must consider the inhomogene-ity of neutron irradiation damage (especially relevant at low fluences) and the formation of B n I m at high N eff,0 through higher order reactions between defects and dopants (at least B i + B ↔ B 2 I). These mechanisms are responsible for the sublinear dependence of φ A on N eff,0, since they affect the efficiency of the removal process and therefore the amount of I mob required to deactivate the dopants. We have analyzed the role of O and C on acceptor removal by including for the formation of VO, C i and B i O defects in our model. The presence of O and C modifies B deactivation mechanisms due to the key role of B i O defects and the trapping of intrinsic defects by O and C. However, for the impurity concentrations used in this work ([O] >> [C]), it has little effect on the overall amount of removed acceptors.
We provide a comprehensive description of N eff for the overall neutron fluences commonly used in high energy physics, by identifying the mechanisms of B deactivation and the simultaneous generation of defect clusters. The analysis at an atomistic level of defect cluster population, combined with the energy levels provided by ab initio and the insight of DLTS/TSC characterization, leads us to propose that the FFC V 4 significantly contributes to N eff at high neutron irradiation fluences.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.