**3. Model implementation**

## **3.1. C++ implementation**

To implement the RWDD model described previously, we choose a C++ programming environment, an object-oriented programming language that offers considerable advantages in terms of advanced structures, such as classes, objects and containers. The C++ code imple‐ menting the RWDD approach was named RWDDCPP.

In this code, the particle track is defined as a C++ container including all the charge packets described as independent objects. The members of the charge packet class include the geo‐ metrical coordinates of the packet and the electrical charge amount per packet. The container content may change during the simulation due to two different mechanisms: (1) minority carrier recombination or (2) carrier extraction that correspond to particles that escape the simulation domain. For modeling the carrier recombination mechanisms, we use a simple exponential law that adjust the number of charges as a function of time: N(t) = N0 exp(−t/τ). In this exponential law, N0 is the initial number of charges deposited by the particle at t = 0 s and the time constant τ is equal to the carrier lifetime.

In RWDDCPP code, the circuit geometry is implemented as a series of 3D rectangular boxes that represent respectively the substrate, the source and drain contacts at the silicon level and the different wells; all these elements constitute the Front-End-Of-Line (FEOL) structure. A simplified Back-End-Of-Line (BEOL) structure may be also modeled, as a stack of insulating material and metal layers. The fine modeling of the reversely biased drain contacts collecting the minority carriers created along the particle track is the most important improvement of the RWDD model at the circuit level. A special "drain class" has been developed, embedding in the same C++ object both the drain contact geometry and the 3D distribution of the electric field induced by the drain bias. At this level, different doping profiles can be taken into account for the p-n junctions: abrupt, gradual or user-defined profiles.

considered individually and all the interesting effects due to particular carriers are summed [16]. In our work, we implemented the first formalism (transport equations) by applying the continuity equation at the collecting (drain) contact. The transient current at the collecting node is directly computed from the number of carriers Δn that reach this contact during the time

In this expression, the displacement current is neglected; that is a reasonable approximation in this case [17]. This collection current is then injected in the electrical simulation to model

Although the model embeds a microscopic physical description of charge generation, transport and collection, the present approach suffers from two evident limitations that should be amended in a future enhanced code version: (1) the calculation is not self-consistent with the electrostatic potential (i.e. the electric field) and (2) it does not implement electrostatic inter‐ actions taking place between charge packets. For this reason, the RWDD model should especially fail in high carrier injection conditions, i.e. at high LETs. These drawbacks will be carefully estimated and quantified in a future work. Another interesting point should be to compare the two methods for computing the collection current to also evaluate the impact of

To implement the RWDD model described previously, we choose a C++ programming environment, an object-oriented programming language that offers considerable advantages in terms of advanced structures, such as classes, objects and containers. The C++ code imple‐

In this code, the particle track is defined as a C++ container including all the charge packets described as independent objects. The members of the charge packet class include the geo‐ metrical coordinates of the packet and the electrical charge amount per packet. The container content may change during the simulation due to two different mechanisms: (1) minority carrier recombination or (2) carrier extraction that correspond to particles that escape the simulation domain. For modeling the carrier recombination mechanisms, we use a simple exponential law that adjust the number of charges as a function of time: N(t) = N0 exp(−t/τ). In this exponential law, N0 is the initial number of charges deposited by the particle at t = 0 s and

Δn I q Δt = ´ (5)

step Δt, i.e.

120 Modeling and Simulation in Engineering Sciences

the circuit response.

**2.4. Model limitations**

the displacement current on transient characteristics.

menting the RWDD approach was named RWDDCPP.

the time constant τ is equal to the carrier lifetime.

**3. Model implementation**

**3.1. C++ implementation**

**Figure 2.** Schematic illustrations of a reversely biased drain junction taken into account in the RWDD approach. The pn junction is defined by its geometry and the bias VR applied to collect and extract carriers from the bulk. A 3D analyti‐ cal model is used to model the space charge region (SCR) in the case of a gradual or an abrupt doping profile. Both the electric field developed in the SCR (E) and the SCR width (W) are controlled by the bias VR. VR (and consequently E and W) may vary in time due to external circuit feedback. (Reprinted with permission from Autran et al. [10], © 2014, Elsevier).

An object of the class "drain" is shown in **Figure 2**; this class is defined by its doping profile (an abrupt junction in this case), its geometrical dimensions (W, L, h), the doping levels of the p and n+ regions (NA and ND, respectively) and the bias (VR) applied to the collecting contact. A 3D analytical model is used to compute the space charge region, since an abrupt p-n junction is considered in this example. The electric field (E) developed in this space charge region and its width (WSCR) are functions of the drain bias potential and can be calculated using the following relations:

$$\mathbf{W}\_{\text{SCR}} = \sqrt{\frac{2\epsilon\_{\text{Si}}\left(\mathbf{V}\_{\text{bi}} + \mathbf{V}\_{\text{R}}\right)}{\mathbf{q}\mathbf{N}\_{\text{A}}}} \tag{6}$$

$$\mathbf{E} = \frac{\mathbf{2}\left(\mathbf{V\_{bi}} + \mathbf{V\_{g}}\right)}{\mathbf{W\_{SC}}} \tag{7}$$

where Vbi is the internal potential of the junction and εSi is the permittivity of silicon.

It must be noted that both WSCR and E can dynamically vary during the transient simulation as a function of the effect of carrier collection on the circuit node potential VR, as illustrated in the following.

#### **3.2. GPU parallelization**

As explained in the previous section, in the RWDD approach, the behavior of each packet of charge is computed independently of the other charges. These processes being independent, the calculation task of the charge transport can be easily parallelized on a graphic processing unit (GPU) whose internal architecture is perfectly adapted to such a massive parallelism. Moreover, RWDD model needs to implement a random number generation procedure that is usually very time-consuming; if the random numbers are independent, this task can also be easily parallelized on GPU. These two characteristics of the RWDD model are expected to offer a considerable computation speed, since the number of parallelizable tasks in RWDD is relatively high. In this work, we used the CUDA programming framework proposed by NVDIA (CUDA version 5.5 64 b) [18]. Random numbers are generated using the CUDA cuRAND random number generation library (Mersenne Twister) within the parallel kernels running on the GPU. **Figure 3** illustrates this parallel implementation of the RWDDCPP code. For the purpose of this work, we separately tested the RWDD algorithm implemented in series on a CPU and in parallel on a GPU. Tests have been performed on a machine equipped with a 3.2 GHz Intel Core i5 and with a NVDIA GeForce GTX 675 MX (960 threads at 600 MHz). Our results show that for the simulation of 100,000 charge packets with 1000 time steps the speedup on GPU was 140× with respect to the classical CPU implementation.

**Figure 3.** Program flow for the parallel version of the RWDDCPP code. The main part is executed on the CPU (host), and the RWDD computational part is executed on the GPU (device). The code is also dynamically coupled with SPICE for circuit simulation within the transient event time domain. (Reprinted with permission from Glorieux et al. [12], © 2014, IEEE).

## **3.3. Multiple node charge collection**

where Vbi is the internal potential of the junction and εSi is the permittivity of silicon.

the following.

2014, IEEE).

**3.2. GPU parallelization**

122 Modeling and Simulation in Engineering Sciences

It must be noted that both WSCR and E can dynamically vary during the transient simulation as a function of the effect of carrier collection on the circuit node potential VR, as illustrated in

As explained in the previous section, in the RWDD approach, the behavior of each packet of charge is computed independently of the other charges. These processes being independent, the calculation task of the charge transport can be easily parallelized on a graphic processing unit (GPU) whose internal architecture is perfectly adapted to such a massive parallelism. Moreover, RWDD model needs to implement a random number generation procedure that is usually very time-consuming; if the random numbers are independent, this task can also be easily parallelized on GPU. These two characteristics of the RWDD model are expected to offer a considerable computation speed, since the number of parallelizable tasks in RWDD is relatively high. In this work, we used the CUDA programming framework proposed by NVDIA (CUDA version 5.5 64 b) [18]. Random numbers are generated using the CUDA cuRAND random number generation library (Mersenne Twister) within the parallel kernels running on the GPU. **Figure 3** illustrates this parallel implementation of the RWDDCPP code. For the purpose of this work, we separately tested the RWDD algorithm implemented in series on a CPU and in parallel on a GPU. Tests have been performed on a machine equipped with a 3.2 GHz Intel Core i5 and with a NVDIA GeForce GTX 675 MX (960 threads at 600 MHz). Our results show that for the simulation of 100,000 charge packets with 1000 time steps the

speedup on GPU was 140× with respect to the classical CPU implementation.

**Figure 3.** Program flow for the parallel version of the RWDDCPP code. The main part is executed on the CPU (host), and the RWDD computational part is executed on the GPU (device). The code is also dynamically coupled with SPICE for circuit simulation within the transient event time domain. (Reprinted with permission from Glorieux et al. [12], © Since the RWDD model has been implemented in an object-oriented language (C++) using dynamic containers, as previously explained in 3.1, complex circuits with an arbitrary number of sensitive areas and collecting nodes can be simulated, intrinsically considering, in this case, multiple node charge collection in the simulation process. The simulation begins with the initialization of the structure and the definition of the circuit nodes corresponding to the different simulated drains. In a second step, the electrical potential of each node is extracted from a steady-state circuit simulation. The values of these potentials are used to initialize the electrical field in the complete structure. Starting the time-domain analysis, at each time step of the transient simulation, the magnitudes of the current sources corresponding to the different collection processes (simultaneously occurring at different circuit nodes) are updated as a function of the number of collected charge packets. **Figure 4** illustrates this multiple node charge collection process in the simple case of two adjacent drain junctions. The time evolution

**Figure 4.** Cartoon (at five different times after the ionizing particle impact) illustrating the time evolution of the charge packets induced by a horizontal ionizing particle impacting two adjacent drains. The collection efficiency is different for these two drains since their biasing state is different (Reprinted with permission from Glorieux et al. [12], © 2014, IEEE.).

of the charge packets induced by a horizontal ionizing particle in this two-node structure is shown at five different times after the ionizing particle impact. As evidenced in **Figure 4**, the collection efficiency of the two drains differs because they correspond to two different nodes in the circuit and their electrical potential (i.e., internal electric field) is not the same in this example.

#### **4. Circuit solving**

Once the current transients have been evaluated using the RWDD model on the different collection nodes of interest in a given device or circuit, one must evaluate the transient electrical response of the circuit. This latter can be directly computed using an internal routine (solving the Kirchoff's circuit laws) or with an external circuit simulator program. We examine here these two solutions, both implemented in our code RWDDCPP.

#### **4.1. Internal routine**

For circuit architectures composed of only a few connected devices, such as CMOS inverters or SRAM cells, the steady-state circuit electrical solution can be simply solved considering Kirchoff's circuit laws and compact models for transistors. In the following, the method is illustrated for the solving of an inverter and a SRAM cell. **Figure 5** (left) defines the different terminal voltages and the source-to-drain current for the NMOS and PMOS transistors; **Figure 5** (right) shows the circuit schematic of a SRAM cell consisting of two cross-coupled CMOS inverters. Considering that the particle strikes the OFF-state NMOS transistor termed NMOS1 (initial conditions V1 = 0, V2 = VDD) the time variations of potential V2 for the sole and isolated inverter #1 can be written as:

$$\frac{\text{d}\mathbf{V}\_{2}}{\text{d}\mathbf{t}} = \frac{-\mathbf{I}\left(\mathbf{t}\right) + \mathbf{I}\_{\text{DVI}}\left(\mathbf{V}\_{1} - \mathbf{V}\_{\text{DD}'}\;\mathbf{V}\_{2} - \mathbf{V}\_{\text{DD}}\right) - \mathbf{I}\_{\text{DNI}}\left(\mathbf{V}\_{1}\;\mathbf{V}\_{2}\right)}{\mathbf{C}\_{\text{N}}} \tag{8}$$

**Figure 5.** Definition of the terminal voltages and the source-to-drain current in n-channel (NMOS) and p-channel (PMOS) transistors (left). Circuit schematic of a SRAM cell formed by two cross-coupled CMOS inverters (right). (Re‐ printed with permission from Autran et al. [10], © 2014, Elsevier).

where I(t) is the current pulse given by Eq. (5) due to the particle strike and collected on the node 2, IDN1 and IDP1 are the currents of the NMOS1 and PMOS1 transistors, respectively, and CN is the node capacitance.

Considering this time, the full SRAM cell composed of the two cross-coupled inverters, we obtain the following system of two coupled differential equations:

$$\begin{cases} \begin{aligned} \frac{\mathbf{dV}\_{1}}{\mathbf{dt}} = & \frac{\mathbf{I}\_{\mathrm{Dr}2}\left(\mathbf{V}\_{2} - \mathbf{V}\_{\mathrm{DD}}, \mathbf{V}\_{1} - \mathbf{V}\_{\mathrm{DD}}\right) - \mathbf{I}\_{\mathrm{DN2}}\left(\mathbf{V}\_{2}, \mathbf{V}\_{1}\right)}{\mathbf{C}\_{\mathrm{N}}} = \mathbf{F}\left(\mathbf{t}, \mathbf{V}\_{1}, \mathbf{V}\_{2}\right) \\\\ \frac{\mathbf{dV}\_{2}}{\mathbf{dt}} = & \frac{-\mathbf{I}\left(\mathbf{t}\right) + \mathbf{I}\_{\mathrm{Dr}1}\left(\mathbf{V}\_{1} - \mathbf{V}\_{\mathrm{DD}}, \mathbf{V}\_{2} - \mathbf{V}\_{\mathrm{DD}}\right) - \mathbf{I}\_{\mathrm{DN1}}\left(\mathbf{V}\_{1}, \mathbf{V}\_{2}\right)}{\mathbf{C}\_{\mathrm{N}}} = \mathbf{G}\left(\mathbf{t}, \mathbf{V}\_{1}, \mathbf{V}\_{2}\right) \end{aligned} \tag{9}$$

Eq. (8) for the sole and isolated inverter #1 or Eqs. (8) and (9) for the full SRAM cell can be easily solved in the time domain, using a fourth-order Runge-Kutta method [19] with a time step Δt identical to the one used in the RWDD algorithm for charge transport and collection. Using notations introduced in Eq. (8), the incremental results for tn+1 = tn + Δt are:

$$\mathbf{V}\_1^{n+1} = \mathbf{V}\_1^n + \frac{\mathbf{K}\_1 + 2\mathbf{K}\_2 + 2\mathbf{K}\_3 + \mathbf{K}\_4}{6} \tag{10}$$

$$\mathbf{V}\_{2}^{n+1} = \mathbf{V}\_{2}^{n} + \frac{\mathbf{L}\_{1} + 2\mathbf{L}\_{2} + 2\mathbf{L}\_{3} + \mathbf{L}\_{4}}{6} \tag{11}$$

with:

of the charge packets induced by a horizontal ionizing particle in this two-node structure is shown at five different times after the ionizing particle impact. As evidenced in **Figure 4**, the collection efficiency of the two drains differs because they correspond to two different nodes in the circuit and their electrical potential (i.e., internal electric field) is not the same in this

Once the current transients have been evaluated using the RWDD model on the different collection nodes of interest in a given device or circuit, one must evaluate the transient electrical response of the circuit. This latter can be directly computed using an internal routine (solving the Kirchoff's circuit laws) or with an external circuit simulator program. We examine here

For circuit architectures composed of only a few connected devices, such as CMOS inverters or SRAM cells, the steady-state circuit electrical solution can be simply solved considering Kirchoff's circuit laws and compact models for transistors. In the following, the method is illustrated for the solving of an inverter and a SRAM cell. **Figure 5** (left) defines the different terminal voltages and the source-to-drain current for the NMOS and PMOS transistors; **Figure 5** (right) shows the circuit schematic of a SRAM cell consisting of two cross-coupled CMOS inverters. Considering that the particle strikes the OFF-state NMOS transistor termed NMOS1 (initial conditions V1 = 0, V2 = VDD) the time variations of potential V2 for the sole and

( ) DP1 1 DD 2 DD DN1 1 2 ( ) ( ) <sup>2</sup>

dV I t I V V , V V – I V ,V

N

**Figure 5.** Definition of the terminal voltages and the source-to-drain current in n-channel (NMOS) and p-channel (PMOS) transistors (left). Circuit schematic of a SRAM cell formed by two cross-coupled CMOS inverters (right). (Re‐


these two solutions, both implemented in our code RWDDCPP.

dt C

printed with permission from Autran et al. [10], © 2014, Elsevier).

example.

**4. Circuit solving**

124 Modeling and Simulation in Engineering Sciences

**4.1. Internal routine**

isolated inverter #1 can be written as:

$$\mathbf{K}\_1 = \mathbf{F}(\mathbf{t}^\mathbf{n}, \mathbf{V}\_1^\mathbf{n}, \mathbf{V}\_2^\mathbf{n}) \,\Delta \,\mathbf{t} \tag{12}$$

$$\mathbf{L}\_1 = \mathbf{G}(\mathbf{t}^n, \mathbf{V}\_1^n, \mathbf{V}\_2^n)\mathbf{\hat{n}}\tag{13}$$

$$\mathbf{K}\_2 = \mathbf{F}(\mathbf{t}^n + \Delta \mathbf{t}/2, \mathbf{V}\_1^n + K\_1/2, \mathbf{V}\_2^n + L\_1/2)\Delta \mathbf{t} \tag{14}$$

$$\mathbf{L}\_2 = \mathbf{G} \{ \mathbf{t}^n + \Delta \mathbf{t} / 2, \mathbf{V}\_1^n + K\_1 / 2, \mathbf{V}\_2^n + L\_1 / 2 \} \Delta \mathbf{t} \tag{15}$$

$$\mathbf{K}\_3 = \mathbf{F} \{ \mathbf{t}^n + \Delta \mathbf{t} / 2, \mathbf{V}\_1^n + K\_2 / 2, \mathbf{V}\_2^n + L\_2 / 2 \} \Delta \mathbf{t} \tag{16}$$

$$\mathbf{L}\_3 = \mathbf{G} \{ \mathbf{t}^n + \Delta \mathbf{t} / 2, \mathbf{V}\_1^n + K\_2 / 2, \mathbf{V}\_2^n + L\_2 / 2 \} \Delta \mathbf{t} \tag{17}$$

$$\mathbf{K}\_{\mathsf{A}} = \mathbf{F}(\mathbf{t}^{n} + \Delta \mathbf{t}, \mathbf{V}\_{1}^{n} + K\_{3}, \mathbf{V}\_{2}^{n} + L\_{3})\Delta \mathbf{t} \tag{18}$$

$$L\_4 = \mathcal{G} (\mathfrak{t}^n + \Delta \mathfrak{t}, \mathcal{V}\_1^n + K\_3, \mathcal{V}\_2^n + L\_3) \Delta \mathfrak{t} \tag{19}$$

Eqs. (10)–(19) constitute the core model implemented in the code RWDDCPP for CMOS inverter and SRAM cell solving. In Eq. (9), NMOS and PMOS source-to-drain currents are analytically modeled using the EPFL-EKV model 2.6 [20, 21], here implemented in C++ as numerical functions (IDN and IDP) having two arguments (see **Figure 5**). The EKV 2.6 MOSFET model is a predictive (scalable) compact model for the simulation of submicron CMOS technologies. It was built taking into account fundamental physical characteristics of the MOS structure. The model offers a continuous modeling of the different regimes of the transistor operation (weak, moderate and strong inversion), which is mandatory for circuit modeling. The version 2.6 takes into account numerous essential issues for the transistor modeling such as the effects of the doping profile, substrate effects, process-related aspects (oxide thickness, effective channel length and width and junction depth), mobility effects due to vertical and lateral electric fields, the velocity saturation, short-channel effects, etc. Section 5 details several examples of simulations using this code RWDDCPP.

#### **4.2. SPICE simulator**

For more complex circuit architectures than a single inverter or a SRAM cell, an external SPICE circuit simulator can be otherwise used in the place of the internal subroutine for circuit solving. In this case, the circuit simulator is instantiated in interactive mode by the executable code with a circuit netlist corresponding to the simulated structure. Current source(s) emu‐ lating the transient collected current(s) at the different circuit nodes is (are) automatically added to the netlist by the RWDDCPP program. At each time step of the simulation, the magnitude of each current source is updated as a function of the number of charge packets collected by the corresponding biased junction, following Eq. (5). Then, the SPICE simulator computes the new voltages on the circuit nodes and finally, from Eqs. (6) and (7), the width of the space charge region(s) and the distribution of the electrical field can be updated, time step by time step. The whole circuit counteraction to the radiation transient can then be taken into account. This approach has been successfully implemented in this work using both NGSPICE [22] and ELDO [23] simulators.
