SITARe: a fast simulation tool for the analysis of disruptive effects on electronics

. This paper is devoted to an exhaustive presentation of a fast computation numerical tool, dedicated to the simulation of transient currents induced by stochastic events in microelectronic devices. This is a part of a numerical platform, SITARe, combining a spice simulator with the semi-analytical model presented here. The paper describes the theoretical model, the calibration. An instance of application illustrates the ability of the tool.


Introduction
With dimensions continuously shrinking, some disruptive effects, as those due to alpha intrinsic impurities, become critical for circuit performances. By testing a circuit in a sheltered underground environment as the LSBB (Laboratoire Souterrain à Bas-Bruit) of Rustrel (France), the alpha contamination level induced by radioactive impurities can be quantified easily. This is a considerable gain of time and precision for the evaluation of an alpha contamination. However, the simulation can provide an interesting complement to this experimental approach. Then the reliability of a circuit can be estimated by simulating a physical disturbance at the scale of an elementary device and by evaluating its consequences at the scale of several interconnected elementary devices.
The simultaneous simulation of the disruptive phenomenon on the elementary device and the electrical behavior of the whole system is hard to perform with the same numerical tool. Then two types of numerical tools are used to simulate the crossing of an ionizing particle in a device, each of them corresponding to a scale of the circuit The first one, the TCAD simulation (as Sentaurus [1]), relies on a quantitative multi physics approach, at the scale of the lowest length in the electronic device. The particle is generated as an electron-hole distribution in the device simulated at a physical level and the transient current is generated at the electrodes by solving the transport equation in the semi-conductor volume [2,3,4]. If well calibrated, TCAD simulation is a realistic method for single events studies but it is time consuming, and generally needs high CPU resources because based on a multi-physics finite elements approach. Finally, TCAD simulation allows a better understanding of the physics at the scale of the elementary device, but is not suitable for an extensive study at the cell scale. The second type of tools relies on the electrical behavior of the entire system or circuit. This approach is based on very classical Kirchhoff's method [5]. Each elementary device is modeled by its electrical biased current characteristic. The impact of the collected currents in one or several nodes is simulated by adding some current sources in the system simulated at circuit level [6,7]. In order to explore all the possibilities in term of time duration or intensity, the currents are randomly generated. The range of value of transient characteristics (duration and level of intensity) of the current can be based on numerical experiments (given by the TCAD simulation for instance) or by the experimental knowledge. The currents can be generated by a simple analytical formulation (a double exponential is often used) and so the physics of disruptive ionizing effect responsible for the current is lost. The results of such simulations allow discussions about the stability or instability of the device under current stresses, giving generally criteria for stability in terms of order of magnitude of the currents and their life time. However, they cannot help in the understanding of the physical links between the nature and the configuration of the disruptive effect and the stability of the cell.
In this context, this paper is devoted to an exhaustive presentation of a fast computation numerical tool allowing the simulation of disruptive effects on electronic devices, from the particle (i.e. alpha) to the loss of functionality. SITARe (Simulation Tool with semi-Analytical method for Transient Effects) connects the simulation of the transient disruptive effect induced by stochastic events, at the scale of a single junction (semi-analytical model), to the simulation at the scale of the electronic device (spice simulator). The disruptive parasitic currents are caused, by an arbitrary localized (in time and in space) generation of electrons-holes pairs near the functional junctions of the device. A semi-analytical model, allowing the fast random generation of statistical data (Monte-Carlo), without losing the physics of the ionizing phenomenon is proposed. The paper describes the theoretical model (section 2), the calibration (section 3) and an instance of application (section 4).

Geometrical set up
The general geometry that can be taken into account by the model is depicted in figure 1. Electron-hole pairs are generated in a volume (called Γ in figure 1) of the silicon bulk. These electric charges move through the volume of the substrate and can be collected at the sensitive surfaces. These surfaces are the junctions of the device (drains or sources for a transistor). Only three are represented in figure 1 (denoted by Σ 1 , Σ 2 , Σ 3 ) as an instance. The number as the shape of collecting surfaces is arbitrary. The tool makes it possible the addition of as collecting surfaces as necessary. Our tool simulates the disruptive currents induced by the creation of electron-hole pairs but not simulate the physical phenomenon at the origin of this generation (e.g. nuclear interaction). The value of ν S is derived from the LET (Linear

Volume charge transport modeling
The volume concentration of the species s (s = h for holes and s = e for electrons) is denoted by ns. ns(M, t) and its associated charge current density vector − → j s (M, t) respect conservation law at any point M in the semi-conductor volume at any instant t. It leads to the local classical general conservation equation: where σ(M, t) denotes the "source term".

Generation and recombination modeling
The right hand σ term of the equation (1) can be considered as a distribution. This approach is classic in electrodynamics, for instance, [6]. It contains the possible recombinations of holes and electrons at M at instant t as the geometry (volume Γ) and the duration of the initial generation of the charges. Let us explain the construction of the term σ we take into account in our model The elementary distribution modeling an instantaneous generation at a point S is modeled as a product between the Dirac distribution peak δ(t = 0) centered at the origin of time (t = 0) and the 3D spatial Dirac peak δ(S ) = δ(x = xS )δ(y = yS )δ(z = zS ) at the point S with coordinates (xS , yS , zS ). The general solution of an arbitrary shape generation in time and space is then obtained by the double convolution between solution obtained for the double Dirac distribution (called Green function) and the temporal shape of the source and the spatial shape of the volume Γ. We choose to model the local recombination of pairs of electron-hole by two distinct local laws of different time range. The first one is a local classical exponential recombination law. It models a recombination rate proportional to the density ns. The second one allows modeling the tail of the current (at long time range) as a power law. It is assumed to be proportional to ns/t. The origin and the validity of this term is discussed after, in a dedicated section. Finally, the σ term is supposed to be as : where νs is a frequency and ζs is a dimensionless coefficient.
Energy Transfer) value. This is the energy lost by the particle, by unit of length, when interacting with the bulk. It varies along the track and is linked to the initial energy of the particle [8]. The values have been extracted from the SRIM tables (Stopping in Range of Ions in Matter) [9]. Interacting with the semi-conductor, this particle deposits electron-holes pairs all along its trajectory with a linear density of ν S per unit of time (S denotes the current point of the track).

Volume charge transport modeling
The volume concentration of the species s (s=h for holes and s=e for electrons) is denoted by n S . n S (M,t) and its associated charge current density vector ȷ " (M,t) respects the conservation law at any point M in the semi-conductor volume at any instant t. It leads to the local classical general conservation equation: where σ (M,t) denotes the "source term".

Generation and recombination modeling
The σ term of the equation (1) can be considered as a distribution. This approach is classic in electrodynamics, for instance [5]. It contains the possible recombination of holes and electrons at M, at the time t as the geometry (volume Γ) and the duration of the initial generation of the charges. The elementary distribution, modeling an instantaneous generation at a point S, is modeled as a product between the Dirac distribution peak δ(t=0) centered at the origin of time t=0 and the 3D spatial Dirac peak δ(S)= δ(x=x s ).δ(y=y s ).δ(z=z s ) at the point S with coordinates (x s , y s , z s ). The general solution of an arbitrary shape generation in time and space is then obtained by the double convolution between the solution obtained for the double Dirac distribution (called Green function) and the temporal shape of the source and the spatial shape of the volume Γ. We choose to model the local recombination of electronhole pairs by two distinct local laws of different time range. The first one is a local classical exponential recombination law. It models a recombination rate proportional to the density n S . The second one allows to model the tail of the current (at long time range) as a power law. It is assumed to be proportional to n S /t. The origin and the validity of this term is discussed after, in a dedicated section. % is a parameter for numerical adjustment. Its sign is a priori unknown, reason we choose to write "+ % ". Finally, the σ term is supposed to be as: where ν S is a frequency and ζ S is a dimensionless coefficient.

Drift diffusion model: link between current and concentration
% , is assumed to be linked to n S by a linear relation combining two process of diffusion in one hand and electrical drift in other hand, at any point M and any time t as: where q S is the electric charge of the species s, D S is the diffusivity of the species s, µ S is the conductivity of the species s, and ξ(M,t) the local electrical field.

Local complete transport equation
Finally introducing relation (3) in (2), it comes a Fokker-Planck equation type [10] as: In the general case, and a priori, µ S as ν S and ζ S could be local and depend on ξ(M,t), M and t. Solving completely the equation (5) in such a context, needs to introduce the conduction bands and then to solve a set of local and non-stationary equations using a spatial mesh of all the structure. The 3D TCAD computation uses this approach. Because our model is based on the will of avoiding such a complexity, some hypothesizes are made to derive an equation with an analytical solution.

Simplifying assumption
1. µ S as ν S and ζ S are supposed to be constant and uniform in the substrate.
2. The local global electric charge is supposed to be null, that is to say that everywhere, the Maxwell-Gauss law writes ∇. , = 0 3. The electric field is not continuous. It is supposed to be null everywhere in the volume of the bulk. At the collecting surface, it is considered as non null, but assumed to be constant in intensity and oriented as the normal of the collecting surface. That is to say, the effect of the space charge zone of the implant is taken into account as if it took place at the implant surface.
With the first hypothesis, the second term of relation (5) writes: The second hypothesis, included in relation (7), leads that in every point M in the volume of the bulk: The third hypothesis, leads that (8), writes, for every point M in the volume of the bulk : Finally these three hypotheses allow to assume that n S (M,t) is driven by the following transport: where D S is the diffusivity of the species s, δ(S), δ(t) are the Dirac distributions centered at the point S and at time t=0, ν S is the density of charges disappearing per unit of time by recombination, supposed constant, the latest term is a recombination term allowing to simulate the tail of the current pulse as a power law (ζ S dimensionless).

Resolution by simplified convolution
The classical solution, of the equation (10), for free space boundaries conditions is called the Green function in free space. It is denoted by g 3d (MS, t) and is given by: where η S is necessarily dimensioned as a frequency.
In the general case, to take into account the spatial shape of the volume Γ, g 3d is convolved with a function describing the spatio-temporal distribution of the generation of electron-hole pairs in the volume Γ. We assume here a last simplifying hypothesis: the process of charges generation is very short compared to the transport time range through the bulk. Then the process is considered as instantaneous and there is no need to operate a convolution on the time variable. ν S0 (M) is the number of element of species s generated per unit of volume around the source point S in the volume Γ. It is dimensioned as the inverse of a volume. Then finally, the solution for n s (M,t) is given by the spatial convolution all along the volume Γ.

Discussion about the recombination modeling
In relation (11), the factor ( % ) S * modulates the classic extinction factor modeled as exp(ν S .t). The global effect of this "law-power" factor is to produce a particular tail (at long time after the maximum of the current is reached) in the transient current shape. It allows to model -with only one coefficient -a part of the complexity of the transport reduced by the simplifying assumptions. Then, the collective transport of charges in a biased medium could be considered as analogous as an anomalous diffusion observed in heterogeneous and complex media. In the diffusion studies, it is usual to consider n s (M,t) as a density of probability of presence and to link the stochastic microscopic rules of displacement (as mobile-Immobile model for instance) to a partial differential equation. More generally, in the stochastic approaches of diffusion problems, the tail (at long time range) is an important characteristic of the associated stochastic process. For instance, this tail is associated to the memory properties of the process. It has been theoretically well explored from many years [11]. It defines the property of stability of the stochastic distribution involved in the transport phenomenon. Stable distributions of stochastic processes have a lot of applications in complex physics problems [12] and signal processing [13]. With the will of avoiding the mathematical complexity of the complete description of stochastic process with memory, we assume the macroscopic factor of modulation in a "law-power" form. The parameter η S is the key of this modulation.

Charge collection
On any surface Σ, described by its local normal vector > Σ, in full generality, the collected current, % V , of the s species, is obtained by integrating the vector % over the surface.
Σ corresponds to all the considered collecting surfaces. In the configuration depicted in figure  1, for instance, it would correspond to the union of Σ 1 and Σ 2 . Because of the relation (4) and the three hypothesis concerning the diffusion current % X (Σ) : where Ω s = μ s .ξ Σ and ξ Σ denotes the scalar amplitude of the electrical field on the collecting surface and normal to the surface.
In Cartesian coordinates, the gradient of n s is easily obtained from relations (11) and (12), that allows estimating the collected current on the surface Σ. One originality of our approach lies in taking into account the two terms in the collected current formulation: diffusion and conduction. The relative modulation of the two terms is ensured by macroscopic independent numerical parameters (Ω s and D s ). For instance, if the collecting surface (the junction) is nonbiased, the diffusion term could ensure alone a non-null current, while Ω s would be chosen as null. Finally, the collected current needs five independent macroscopic parameters to be modelled: D s , Ω s , ν s , η s , ζ s . Their numerical values are fixed at the end of a numerical calibration process. One instance is given in section 3.
The quantity of charge % (Σ, N ), passing through the collecting surface from the beginning of the collection t i to the time t 0 is evaluated as: Finally, the modeling of the junction is numerically calculated in a classical way by discretizing the surface in rectangular blocks. The only originality comes from an adaptive time step, based on a a priori and not a posteriori calculation. It is optimized to get a discretization around the current peak and release the time step in the relaxation part. The time T b corresponding to the maximum of current depends on the minimal distance between the track and the junction, r min . Its obtained from % = 0 as:

Calibrating the computation
The calibration of our code consists in adjusting the parameters of our model to fit as well as possible the intensity given by a reference computation, here provided by a TCAD simulator (Sentaurus Device). It is a 3D Finite Elements Method commercial code taking into account the physics of the semiconductor in a nonstationary mode. All the quantities (n s , etc..) are locally and temporally computed, point by point and time after time. TCAD simulation needs the entire knowledge of the physical parameters in the bulk (doping profiles, junction dimensions, Because the collecting surface (junction) can be externally polarized or not, at least two TCAD simulations (polarized or not) are needed to determine the "fit parameters". Then few (one, two or three) reference-configurations for Γ are used to refine the fitting process. To illustrate this latest, one instance of results of calibration is presented in figure 3. The volume Γ is an horizontal line, located at 0.5 µm underneath the collecting surface, centered on the center of the implant (the middle of the track is in front of the center of the implant), and oriented in the same direction than the smallest dimension of the implant (x axis). It is 0.5 µm long. It corresponds to an impact of an α particle, ionizing the bulk, with a constant LET (linear energy tranfser) corresponding to 1.5 MeV. Those values are extracted from SRIM database [25]. The collecting surface is an N doped junction of 0.32 µm large and 1.1 µm long. This implant can be polarized (configuration C1) or not (configuration C2) ). The nature of the collected charges is also electrons (e for the subscripts). In this example, it is aimed to differentiate the two configurations with only the conduction parameter ( Ωe). That is to say, the calibration process is enforced to have same charges transport physics in the volume for the two configurations (νe and De are the same). Thus, between the two configurations, there are three degrees of freedom (ηe, ζe, Ωe).  A very light shift on position of the maximum of the currents (∼ 2 × 10 −2 ns) is observed in C1 configuration between TCAD and semi analytical results. On the other hand, in the C2 configuration the major difference between the two curves appears in the the tail of the curve, between 0.2 and 0.4 ns. Nevertheless the global agreement seems to be good. The choice of the level of agreement and the areas where the fit has to be optimized is determined by the aimed application. For instance one classical approach is to fit the Imax, tmax area, where the current is maximum. Because the collecting surface (junction) can be externally polarized or not, at least two TCAD simulations (polarized or not) are needed to determine the "fit parameters". Then few (one, two or three) reference-configurations for Γ are used to refine the fitting process. To illustrate this latest, one instance of results of calibration is presented in figure 3. The volume Γ is an horizontal line, located at 0.5 µm underneath the collecting surface, centered on the center of the implant (the middle of the track is in front of the center of the implant), and oriented in the same direction than the smallest dimension of the implant (x axis). It is 0.5 µm long. It corresponds to an impact of an α particle, ionizing the bulk, with a constant LET (linear energy tranfser) corresponding to 1.5 MeV. Those values are extracted from SRIM database [25]. The collecting surface is an N doped junction of 0.32 µm large and 1.1 µm long. This implant can be polarized (configuration C1) or not (configuration C2) ). The nature of the collected charges is also electrons (e for the subscripts). In this example, it is aimed to differentiate the two configurations with only the conduction parameter ( Ωe). That is to say, the calibration process is enforced to have same charges transport physics in the volume for the two configurations (νe and De are the same). Thus, between the two configurations, there are three degrees of freedom (ηe, ζe, Ωe).  A very light shift on position of the maximum of the currents (∼ 2 × 10 −2 ns) is observed in C1 configuration between TCAD and semi analytical results. On the other hand, in the C2 configuration the major difference between the two curves appears in the the tail of the curve, between 0.2 and 0.4 ns. Nevertheless the global agreement seems to be good. The choice of the level of agreement and the areas where the fit has to be optimized is determined by the aimed application. For instance one classical approach is to fit the Imax, tmax area, where the current is maximum. To give 8 gate oxide thickness…), usually given by the manufacturer of the simulated device [14]. Then the layout allows to define roughly the geometrical function z = f (x, y) for the collecting surfaces. The fitting parameters are then the physical quantities of the substrate (D s , ν s , η s , ζ s ), and the global parameter, Ω s (relation 14). Because the collecting surface (junction) can be externally biased or not, at least two TCAD configurations (biased or not) are needed to determine the fitting parameters. Then few (one, two or three) reference-configurations for Γ are used to refine the fitting process. The most convenient cases, for the calibration of diffusion process (track generated out of the junction), are the vertical impact (less favorable diffusion case) and a horizontal impact (most favorable diffusion case). One instance of calibration is presented in figure 2. The volume Γ is a 0.5 μm horizontal line, located at 0.5 μm underneath the collecting surface, centered on the center of the junction, and oriented in the same direction than the smallest dimension of the implant (x axis). It corresponds to an α particle impact, ionizing the bulk, with a constant LET corresponding to an initial energy of 1.5 MeV. The collecting surface is a N doped junction (0.32 μm by 1.1 μm). This junction can be biased (configuration C1) or not (configuration C2). The collected charges are electrons (e for the subscripts). In this example, the two configurations are differentiated only with the conduction parameter (Ω e ). It means that the calibration process is enforced to have same charges transport physics in the volume for the two configurations (ν e and D e are the same). Thus, between the two configurations, there are three degrees of freedom (η e , ζ e , Ω e ). Note that the TCAD calculation time is around 48 hours, in a standard workstation computer while it is less than 1 second for the semi-analytical computation.

One instance of application
The figure 3 presents one instance of application for alpha tracking [15]. 9 junctions (figure 3(a)) are simulated by using SITARe. The 9 corresponding currents extracted from SITARe (figure 3(b)) are injected in 9 voltage controlled oscillators used for particle detection. By treating the 9 signals at the output of the 9 detectors, the crossing of the alpha particle through the 9 by 9 matrix of detectors can be reconstructed. This example illustrates the ability of SITARe to study the physical behavior linked to a a given device. Another example published in [14], illustrates the capacity of the tool to treat a huge number of data in order to extract simplified metrics.

Discussion and conclusion
In this paper, an intermediate numerical tool able to simulate the collection of charges, generated by a radiative particle, on several biased junctions is exposed. Similar approaches have been developed since many years for the study of radiation effects in electronic devices. Usually, these tools are dedicated to Single Event Rates prediction in a given radiative environment [7,8]. In comparison, our tool does not aim at prediction but at the analysis of a huge number of data to extract significant metrics for reliability prediction and the analysis of the physical behavior of any device.
The main limitation of our tool lies in the calibration process. We have to launch TCAD 3D complete and transient simulations. This multi physics simulation needs a lot of information directly provided by the manufacturer of the technology. These information is often confidential and not easy to obtain. At this time, this point is difficult to overpass because whatever the involved model, or its numerical implementation, the calibration process needs to have electric characteristics used as reference such as described in the previous section.
Another limitation lies in the huge quantity of data generated for a statistical study. In SITARe, this quantity is reduced before the electrical simulation. A post treatment is necessary after the electrical simulation to keep results easy to exploit. For a single collecting junction, one post-treatment has been successfully proposed in [16]. However, when several collection surfaces are involved, the solution proposed becomes ineffective because there are as many diagrams as collecting surface. We are working on a systematic reduction of the representation of the experiences, before and after electrical simulations, by considering some pertinent physical quantities such as energy or action at the device scale, and not only at the junction scale.