Compact SER Models via Model Order Reduction of Diffusion-Based Charge Collection

Compact SER Models via Model Order Reduction of Diffusion-Based Charge Collection Pavlos Stoikos, Olympia Axelou, Pelopidas Tsoumanis, Georgios-Ioannis Paliaroutis, and George Floros Department of Electrical and Computer Engineering, University of Thessaly, Volos, Greece {pastoikos, oaxelou, petsouma, gepaliar, gefloros}@e-ce.uth.gr Abstract—In contemporary VLSI designs, reliable Soft Error Rate (SER) estimation is critical to ensure circuit reliability under ionizing radiation. Physics-based SER models leverage the inherent linearity of the diffusion equation to construct a discrete representation using the Finite Difference Method (FDM) in three dimensions. However, traditional methods often rely on detailed numerical simulations of these physical processes, resulting in high computational costs and limited scalability across various technological nodes. To accelerate simulation process, in this paper, we propose a compact SER modeling approach based on Moment Matching (MM) combined with the Extended Krylov Subspace (EKS). Experimental results demonstrate significant computational speedups, up to 21×, while maintaining the accuracy of the full model. Index Terms—Soft Error Rate, Diffusion-Collection Equation, Model Order Reduction, Krylov Subspace I. INTRODUCTION With the ongoing CMOS technology scaling down to advanced nanometer nodes, circuits in critical applications like aerospace, automotive, and data centers become increasingly susceptible to soft errors, i.e., transient faults caused by radiation-induced charge deposition. Soft errors are primarily triggered by high-energy particles, such as cosmic rays, alpha particles, and atmospheric neutrons, which generate electronhole pairs when they strike a semiconductor substrate [1]. Furthermore, as supply voltages decrease and node capacitances shrink, the critical charge, Qcrit, required to trigger an error also reduces, making Soft Error Rate (SER) an increasing concern in both terrestrial and space applications. As a result, accurate modeling and efficient prediction of SER are essential for designing radiation-hardened circuits and ensuring the reliability of CMOS-based systems [2]. A high-energy particle striking a semiconductor device initiates a series of physical phenomena that can potentially result in a soft error [3]. As the incident particle traverses the material, it gradually loses energy, disrupting the semiconductor lattice, and creating a dense cylindrical trail of electron-hole pairs. Subsequently, these charge carriers are rapidly transported and collected by sensitive nodes via drift and diffusion mechanisms, driven by inherent electric fields. In a memory cell, the resultant transient current may cause a bit-flip, whereas in a logic cell, such disturbance manifests as a glitch at the gate output, known as a Single Event Transient (SET). An SET can propagate and induce a soft error if it is latched by a sequential element [2]. To capture the *This work has been done in the framework of EU-funded Horizon Europe Twinning project TWIN-RELECT, under the Grant Agreement No. 101160314. impact of ionizing radiation on IC operation, it is essential to accurately model these underlying physical processes, putting emphasis on solving the diffusion equations, thus, enabling and facilitating cross-layer reliability assessments from the device to the system level. While the full analytical approach based on matrix exponential provides an exact solution to the charge transport equations, it becomes computationally expensive for large-scale semiconductor systems due to the high dimensionality of the state-space [4]. Similarly, the Random-Walk Drift-Diffusion (RWDD) method suffers from excessive computational cost and noise due to the need for tracking a large number of particles over time [5]. As a result, Model Order Reduction (MOR) offers a scalable alternative by significantly reducing computational complexity while preserving the dominant charge collection dynamics, enabling efficient and accurate SER estimation across different CMOS process nodes. In this paper, we present a novel MOR-based methodology capable of extracting accurate and compact SER models. The main contributions of this work are summarized as follows. Firstly, we apply an optimized Moment Matching (MM) approach using the Extended Krylov Subspace (EKS) projection [6], which significantly reduced the computational complexity while preserving accuracy. Secondly, we incorporate Monte Carlo simulations to account for statistical charge deposition variations, enhancing the reliability of SER predictions. Finally, experimental results demonstrate that our methodology achieves orders-of-magnitude speedup compared to the full analytical model, while maintaining high accuracy. The rest of this paper is organized as follows. Section II provides a brief overview of the physics-based diffusioncollection SER modeling. In Section III, we present our methodology, which is a computationally efficient approach for extracting compact SER models using MM and Monte Carlo simulations. Section IV presents our experimental results. Finally, conclusions are drawn in Section V. II. PHYSICS-BASED SER MODELING Electronics operating at terrestrial and space levels are constantly exposed to cosmic rays and energetic particles that can potentially induce soft errors by disrupting either memory data or circuit proper operation. A soft error is the result of complex physical phenomena occurring as a high-energy particle strikes and interacts with a semiconductor, usually a silicon device, generating excess charge through direct or indirect ionization [7]. Messenger [3], based on the physics behind these phenomena, developed an approximate analytical solution for the generated current, which has been widely used until now. A. Diffusion–Collection Mechanism The number of generated electron-hole pairs δn0per unit length depends on the device characteristics and the ionizing particle’s Linear Energy Transfer (LET), which is the energy per unit length that the particle deposits to the material along the track, and is given by: δn0=Eparticle Ee−h (1) where Eparticle is the energy that the particle deposits in the device and Ee−his the minimum energy required to generate a single electron-hole pair in silicon (i.e., 3.6eV ). The total generated charge Qtotal along a track with length Lis calculated as follows: Qtotal =qδn0L=qEparticle Ee−h L(2) where qis the elementary charge (i.e., 1.6×10−19C). After energy deposition and charge generation, the charge carriers undergo transport governed by drift and diffusion. The fundamental diffusion governing equation is [3]: ∂δn(r, t) ∂t −δn(r, t) τ=D∗∇2δn(r, t)(3) where τis the carrier lifetime and D∗is the ambipolar diffusion coefficient, which at very high carrier concentrations is given by: D∗= 2DnDp/(Dn+Dp)(4) where Dnis the electron diffusion coefficient and Dpis the hole diffusion coefficient. Eq. (3) describes how the excess charge carriers spread out and recombine over time tfor all points within the material being at radial distance rfrom the ion track. B. Discretization of the Diffusion Equation in 3D A common procedure for the numerical solution of Eq. (3) involves discretizing the three spatial coordinates with steps ∆x,∆y, and ∆z. The diffusion equation in three dimensions is initially written as: ∂δn(r, t) ∂t −δn(r, t) τ=D∗(∂2δn(r, t) ∂x2+∂2δn(r, t) ∂y2+∂2δn(r, t) ∂z2) (5) By approximating the second-order derivatives using central finite differences, we derive the following expression for δni,j,k, representing the number of pairs at each discrete point (i, j, k)with its neighboring points: dδni,j,k dt =D∗(δni+1,j,k −2δni,j,k +δni−1,j,k ∆x2 +δni,j+1,k −2δni,j,k +δni,j−1,k ∆y2 +δni,j,k+1 −2δni,j,k +δni,j,k−1 ∆z2)−δni,j,k τ (6) This discretization process results in a system of ordinary differential equations (ODE) given by: ˙ δn(t) = Aδn(t)(7) with the corresponding matrix elements for a single grid point (i, j, k)is defined as follows: A(i,j,k),(i,j,k)=−2D∗ ∆x2+D∗ ∆y2+D∗ ∆z2+1 τ and for its six nearest neighbors: A(i+1,j,k),(i,j,k)=A(i−1,j,k),(i,j,k)=D∗ ∆x2 A(i,j+1,k),(i,j,k)=A(i,j−1,k),(i,j,k)=D∗ ∆y2 A(i,j,k+1),(i,j,k)=A(i,j,k−1),(i,j,k)=D∗ ∆z2 Additionally, in order to incorporate the external charge injection term Bu(t), we define Eq. (7) as the following Linear time-invariant (LTI) system: ˙ δn(t) = Aδn(t) + Bu(t)(8) where Bspecifies the spatial locations where the external source u(t)injects charge, while u(t)represents the external charge generation. The external charge generation function u(t)is given by: u(t) = (δn0, t = 0 0, t > 0(9) C. SER Evaluation The collected current at a sensitive node considering a unique point charge is given by: Icol(r, t) = qASvcolδn(r, t)(10) where area ASis the sensitive drain area and the vcol is the average carrier velocity in the space charge region [4], [8]. By integrating Icol(r, t)over time, the total collected charge is obtained: Qcol(r, t) = Zt 0 Icol(r, t′)dt′(11) However, if we consider discrete time points, the integration can be approximated numerically using a trapezoidal summation: Qcol(r, t) = N X i=1 Icol(ti)·∆ti(12) The primary purpose of this work is to assess the current pulses and the collected charges resulting from particle hits, substituting TCAD simulations. In [9], these parameters have been evaluated by operating TCAD simulations. Even though this process is regarded as highly accurate, on some occasions it can be time-consuming and computationally complex. Thus, in this paper, an alternative approach is proposed to estimate the aforementioned factors, maintaining accuracy and decreasing computational overhead. Following the methodology presented in [9], the next step involves applying the extracted current pulses to gate nodes in SPICE simulations, modeling the SET pulse characteristics, its width and height. The final target is to evaluate the circuit SER by integrating the Monte Carlo method for different primary inputs scenarios, SET propagation analysis, and masking mechanisms modeling into the SER evaluation tool. III. COMPACT SER MODELS A. MM Overview For the LTI system defined in Eq. (8), the aim of MOR is to create a reduced model: ˜ ˙ δn(t) = ˜ A˜ δn(t) + ˜ Bu,˜y(t) = ˜ L˜ δn(t)(13) where ˜ A∈Rr×r,˜ B∈Rr×m, and ˜ L∈Rp×rwith r≪n. Note that Lrepresents the charge collected at the contact points. By applying the Laplace transform to the system and assuming the initial condition ∆n(0) = 0 with an impulse response U(s)=1, the system can be rewritten as: (sI−A)∆n(s) = B,Y(s) = L∆n(s)(14) Considering the Taylor series expansion of X(s)around zero, the transfer function of Eq. (8), i.e., H(s) = L(sI− A)−1B, can be expanded in terms of its moments as follows: H(s) = M0+M1s+M2s2+M3s3+. . . (15) where M0,M1,M2,M3, and so on, represent the moments of the transfer function. The goal of MM techniques is to create a reduced-order model (ROM) where certain moments of the reduced-order transfer function, ˜ H(s), denoted as ˜ Mi, align with corresponding moments from the original transfer function H(s). B. EKS Subspace To generate the reduced-order model, MM methods compute an iterative projection subspace and project the original system onto this subspace for extraction. Specifically, for a reduced model of order r, and the number of moments k=r p, the matrix V∈Rn×r(r≪n) is the projection matrix that spans the k-dimensional Krylov subspace: range(V)=span{B,AB,...,AkB} The reduced-order model is then computed as follows: ˜ A=VTAV,˜ B=VTB,˜ L=LV (16) where ˜ A∈Rr×r,˜ B∈Rr×p, and ˜ L∈Rq×r. The reduction process itself is independent of how the subspace is selected, however, its success mostly depends on the chosen subspace. In order to address this, the EKS subspace can be utilized, particularly because determining the reduced-order model around multiple expansion points lacks a unified method and varies by the different problem conditions. EKS combines the standard Krylov subspace, previously described, with the inverse Krylov subspace derived from A−1. Consequently, the projection subspace can be computed as follows: range(V)=span{B,A−1B,AB,...,A(m/2)−1B,A−m/2B} (17) Finally, the projection matrix V∈Rn×r, which spans the EKS subspace defined in Eq. (17), can be computed using any Arnoldi procedure. Specifically, the calculation starts with [B,A−1B]and iteratively expands the EKS subspace until an ”a-priori” selection of moments is reached. C. Monte Carlo Ionizing radiation leads to charge deposition in semiconductor devices. It is inherently stochastic and follows a statistical distribution, influenced by factors such as particle energy, material properties, and incident position. By modeling this uncertainty, we obtain a robust SER prediction that captures real-world charge transport variability. However, directly solving the diffusion equation for each Monte Carlo iteration is computationally expensive; thus, we leverage MOR techniques to accelerate the simulation while preserving the accuracy of charge transport dynamics. Previous works have demonstrated that modeling charge deposition as a Gaussian distribution effectively characterizes charge collection in semiconductor detectors exposed to ionizing radiation [10]. This approach simplifies Monte Carlo sampling and enables efficient SER estimation while capturing realistic charge deposition variability. The charge deposition per radiation event, Qdep, can be expressed as: Qdep ∼ N(µQ, σQ)(18) where µQrepresents the mean deposited charge and σQ denotes its standard deviation. As a result, in our approach, ion strike is modeled as a Gaussian distribution of charge carriers. IV. EXPERIMENTAL RESULTS To evaluate the efficiency of the proposed compact SER models, we designed a set of silicon-based semiconductor benchmarks. Specifically, we defined the semiconductor properties used in our experimental analysis. We assumed a diffusion coefficient of D∗= 533cm2/s, which corresponds to high-mobility charge carriers in silicon. The charge collection velocity was set to ucol = 3.55 ×106cm/s, modeling the drift of carriers toward a sensitive node under an applied electric field. The charge injection was modeled as a localized point source, and for each benchmark, we used a different contact area and distance from the deposition site. For all experiments, we assumed that the initial deposition was set at (Lx/2,0,0), where Lx/2represents the x-coordinate size in µm, and the injection was initialized with 50,000 electrons in the first two benchmarks and as a Monte Carlo approach in the remaining benchmarks. The proposed approach was implemented in Matlab 2024.1 with the built-in packages for numerical analysis. Our experiments were performed on a Windows workstation with 32 GB memory and a 4.7 GHz Intel Core i9-9900k processor with 16 threads. To demonstrate the scalability and precision of our method, we applied transient analysis on both the original model and the ROM, using the Backward Euler method and simulating up to 45ps. From Table I, we observe that our method achieves TABLE I RUNTIME AND ACCURACY COMPARISON BETWEEN THE ORIGINAL AND REDUCED MODEL FOR LARGE-SCALE 3D BENCHMARKS SER models Size (µm3) Discr. step (µm) Contact area Original model Reduced model Speedup Accuracy A (µm2) Dist. (µm) Order Runtime (s) Order Runtime (s) Bench. #1 1x1x1 0.09 0.02x0.02 1.02 1331 0.008 13 0.0005 15.54×2.30E-03% Bench. #2 3x3x3 0.08 0.5x0.5 3.15 54872 0.901 28 0.046 19.72×4.50E-04% Bench. #3 6x4x4 0.09 0.6x0.4 4.47 129712 4.395 34 0.214 20.56×1.93E-03% Bench. #4 8x8x8 0.1 1.5x1.2 8.52 512000 58.784 52 4.288 13.71×1.61E-02% Bench. #5 10x10x8 0.09 1.5x1.5 3.15 1096569 420.809 83 21.411 19.65×2.41E-02% a maximum runtime speedup of 21×, while maintaining a relative error below 0.02%, even for the most computationally demanding benchmarks. Furthermore, we achieve a system matrix size reduction up to 99.5%, demonstrating the efficiency of our approach. Moreover, Fig. 1 illustrates the currents from Eq. (10) and the number of collected electrons defined by Eq. (12) at the contact area of benchmark#2. Both methods yield almost identical profiles, indicating the proposed method effectively predicts charge collection dynamics over the time. Finally, Fig. 2 presents the carrier density distribution at t= 0.4ps and t= 3ps in the (x, y, 0) plane of benchmark#2. Specifically, the scatter plots depict the charge particles moving due to diffusion, while the colormaps represent the charge density distribution in space. From the plot’s behavior, we conclude that our method accurately captures charge transport according to the diffusion-collection model. Time (ps) Current (Original) Charge (Original) Charge (Reduced) Current (Reduced) 60 40 20 80 100 120 140 0 25 Current (μA) 30 20 10 40 60 70 80 0 20 30 10 5 35 40 Charge (Electrons) Fig. 1. Time evolution of current and collected charge at the contact point for Bench. #2. V. CONCLUSION In this paper, we have presented an efficient methodology for extracting compact SER models, leveraging both MM and Monte Carlo simulations. Starting with the physics-based diffusion-collection SER modeling approach, we have demonstrated how reduced-order models can effectively predict the charge collection dynamics at points of interest. It can be observed by experimental results, that our approach can achieve 18×average runtime speedup with negligible relative error. REFERENCES [1] R. C. Baumann, “Radiation-induced soft errors in advanced semiconductor technologies,” IEEE Transactions on Device and materials reliability, vol. 5, no. 3, pp. 305–316, 2005. 1.0 0.5 0.0 1.5 2.0 2.5 3.0 3.15μm 0.5 1.5 2.5 1.0 0.0 3.0 2.0 t = 0.4 ps (As =0.5×0.5μm2) Contact area Carrier density 1.0 0.5 0.5 0.0 1.5 1.5 2.0 2.5 2.5 3.0 ×103 1.0 0.0 3.0 2.0 5 15 25 35 t = 0.4 ps 1.0 0.5 0.0 1.5 2.0 2.5 3.0 0.5 1.5 2.5 1.0 0.0 3.0 2.0 t = 3 ps Carrier density 1.0 0.5 0.5 0.0 1.5 1.5 2.0 2.5 2.5 3.0 ×10 1.0 0.0 3.0 2.0 20 40 60 80 t = 3 ps Fig. 2. Charge distribution in the projected xy plane (z= 0) of Bench. #2 at 0.4 ps and 3 ps: Scatter plots (left) and carrier density maps (right) show charge diffusion and drift toward the contact area (3.15µm away). [2] G. I. Paliaroutis, P. Tsoumanis, N. Evmorfopoulos, G. Dimitriou, and G. I. Stamoulis, “A placement-aware soft error rate estimation of combinational circuits for multiple transient faults in cmos technology,” in 2018 IEEE International Symposium on Defect and Fault Tolerance in VLSI and Nanotechnology Systems (DFT). IEEE, 2018, pp. 1–6. [3] G. Messenger, “Collection of charge on junction nodes from ion tracks,” IEEE Transactions on nuclear science, vol. 29, no. 6, pp. 2024–2031, 1982. [4] J.-L. Autran and D. Munteanu, “Physics-based analytical formulation of the soft error rate in cmos circuits,” IEEE Transactions on Nuclear Science, vol. 70, no. 5, pp. 782–791, 2023. [5] M. Glorieux, J. L. Autran, D. Munteanu, S. Clerc, G. Gasiot, and P. Roche, “Random-walk drift-diffusion charge-collection model for reverse-biased junctions embedded in circuits,” IEEE Transactions on Nuclear Science, vol. 61, no. 6, pp. 3527–3534, 2014. [6] G. Floros, N. Evmorfopoulos, and G. Stamoulis, “Frequency-limited reduction of regular and singular circuit models via extended krylov subspace method,” IEEE Transactions on Very Large Scale Integration (VLSI) Systems, vol. 28, no. 7, pp. 1610–1620, 2020. [7] P. E. Dodd and L. W. Massengill, “Basic mechanisms and modeling of single-event upset in digital microelectronics,” IEEE Transactions on nuclear Science, vol. 50, no. 3, pp. 583–602, 2003. [8] S. Uznanski, Single Event Upsets in Sub-65nm CMOS technologies. LAP LAMBERT Academic Publishing, 2011. [9] P. Tsoumanis, G. I. Paliaroutis, N. Evmorfopoulos, and G. Stamoulis, “Analysis of the impact of electrical and timing masking on soft error rate estimation in vlsi circuits,” Technologies, vol. 10, no. 1, p. 23, 2022. [10] G. R. Srinivasan, “Modeling the cosmic-ray-induced soft-error rate in integrated circuits: An overview,” IBM Journal of Research and Development, vol. 40, no. 1, pp. 77–89, 1996.