Victor Malherbe

Gilles Gasiot

Thomas Thery

Jean-Luc Autran

Philippe Roche

Accurate Resolution of Time-Dependent and Circuit-Coupled Charge Transport Equations: 1D Case Applied to 28 nm FD-SOI Devices

We present a 1D drift-diffusion solver for singleevent simulation. Owing to its computational speed and circuitcoupling ability, the module is embedded in our soft error rate simulation platform, enabling projections on logic cells in 28 nm FD-SOI.

I. INTRODUCTION

The physics of radiation-induced charge transport and collection is one of the key blocks in the analysis of single event effects on semiconductor devices. Because the constituting equations of electrostatics and carrier motion are tightly coupled, the formulated problem is inherently stiff and best investigated within Technology Computer Aided Design (TCAD) simulations. However, such numerical codes are computationally intensive, typically requiring several hours for a single impact simulation on a 3D structure. Thus TCAD is neither suitable for error rate prediction nor interactive feedback at circuit-design time: both of them require single event models compact enough to execute in seconds, yet physically accurate. Within the fuzzy notion of "physical accuracy", two peculiar features ought to be tackled by such models if they are to be used at analog-circuit ("SPICE") level: time dependence, and circuit coupling. The Ambipolar Diffusion with Cut-off (ADC) model [START_REF] Edmonds | A Theoretical Analysis of Steady-State Charge Collection in Simple Diodes Under High-Injection Conditions[END_REF] offered extremely valuable insight into the collection mechanisms, but the analysis was carried out under steady-state conditions only. That same idea of considering purely ambipolar transport in the absence of electric field, followed by collection from the electric field in a depletion region, gave rise to the "diffusion-collection" model [START_REF] Palau | Device simulation study of the SEU sensitivity of SRAMs to internal ion tracks generated by nuclear reactions[END_REF]: the time-dependent diffusion equation can be solved analytically, however due to the separation between diffusive and drift transport, no true circuit coupling can be implemented. The authors of [START_REF] Kauppila | A Bias-Dependent Single-Event Compact Model Implemented Into BSIM4 and a 90 nm CMOS Process Design Kit[END_REF] were among the first to tackle the question of circuit coupling, replacing the (independent) current source by a bias-dependent current-generating network comprised of a capacitor and a (dependent) current source. Circuit-element formulations have also proved to be popular for SOI technology, where bipolar amplification (i.e. additional injection due to the disturbed body potential) is commonly accounted for by injecting the primary radiation-induced current into the base of a bipolar junction transistor [START_REF] Fulkerson | Prediction of SOI single-event effects using a simple physics-based SPICE model[END_REF]. Equivalent-circuit approaches like these can offer great predictive power but only after proper calibration, which is perhaps their main weakness when compared to more physics-based models. In this paper we develop a single event model based on direct resolution of the semiconductor equations on a onedimensional mesh. The thrust behind this "1D" assumption is twofold: firstly, it greatly alleviates computational cost, allowing to meet the speed requirement for SPICE-level single event models. Secondly, it is physically relevant to advanced technologies such as Fully-Depleted SOI (FD-SOI) or FinFET, which both exhibit active silicon volumes whose dimensionality can be reduced: unlike bulk technology where the transport problem is intrinsically three-dimensional, in 28 nm FD-SOI the "thickness axis" can actually be collapsed thanks to the ultra-thin (7 nm) silicon film, and the "width axis" may also be neglected due to the electric field's orientation from drain to source. Last but not least, coupling with an external circuit is enabled thanks to the interactive mode of Mentor Graphics' Eldo circuit solver [START_REF]Eldo Classic -Foundry Certified SPICE Accurate Circuit Simulation[END_REF], which makes it possible to run a transient simulation step by step with dynamic access to the circuit parameters. The paper is organized as follows: we first go through the equations that need to be solved (including how to discretize them properly), and shed some light on practical implementation. We then validate the model against more detailed TCAD simulations, finally to integrate it inside our Monte-Carlo radiation simulator to perform cell-level soft-error studies in 28 nm FD-SOI.

II. SOLVER FORMULATION A. Physical models

Guided by Occam's razor, we wish to only include in our solver those physical models that are necessary to account for single event mechanisms, especially in FD-SOI. We thus consider the drift-diffusion transport equations (e.g. instead of more sophisticated hydrodynamic transport):

∇ • εF = q(p -n + N d -N a ) (1)
∂n ∂t + R -G - 1 q ∇ • J n = 0 (2) ∂p ∂t + R -G + 1 q ∇ • J p = 0 (3)
J n = qµ n nF n + qD n ∇n (4)
J p = qµ p pF p -qD p ∇p (5)
(1) is Poisson's equation, where ε is the dielectric permittivity of silicon, F = -∇V is the electric field obtained from electrostatic potential V , q is the elementary charge, n (resp. p) is the electron (resp. hole) density, and N d (resp. N a ) is the donor (resp. acceptor) density assumed to be totally ionized. ( 2) and ( 3) are the electron and hole continuity equations, where R and G are the recombination and generation terms, and J n (resp. J p ) is the electron (resp. hole) current density. Finally (4) and ( 5) are the drift-diffusion equations for electrons and holes, where µ n , D n , µ p , D p are their respective mobilities and diffusivities. The driving fields F n,p = -∇V n,p = -∇(V ± 1 2 V bgn ) account for the spatial variations of both the electrostatic potential and the energy gap E g,ef f = E g -qV bgn , the latter arising because of bandgap narrowing (BGN) at high doping densities. We account for BGN with the del Alamo model [START_REF] Del Alamo | Simultaneous measurement of hole lifetime, hole mobility and bandgap narrowing in heavily doped n-type silicon[END_REF], and the 1 2 factors imply equal splitting between the conduction and valence bands. We assume Maxwell-Boltzmann statistics for the carrier concentrations (as opposed to Fermi-Dirac which are more tedious to implement), because our TCAD simulations only showed small differences in terms of radiation currents with or without Fermi statistics (and therefore, the mobilities and diffusivities are related by Einstein's relation: D = µk B T /q). On the other hand, note that overlooking bandgap narrowing leads to completely unrealistic predictions for the bipolar amplification mechanism: when the mass-action law is modified into (np) eq = n 2

int,ef f = n 2 int exp(qV bgn /k B T ) where n int is the intrinsic carrier concentration, then with highly-doped implants the minority concentration is off by an exponential factor, yielding overly pessimistic projections for the parasitic bipolar gain and thus the total collected charge under radiation events. We also model temperature-dependent bandgap and conduction/valence band effective densities of states [START_REF] Bludau | Temperature dependence of the band gap of silicon[END_REF], [START_REF] Green | Intrinsic concentration, effective densities of states, and effective mass in silicon[END_REF], [START_REF] Lang | Temperature dependent density of states effective mass in nonparabolic p-type silicon[END_REF]. Finally, note that we define the electrostatic potential with respect to the intrinsic Fermi level, and therefore at equilibrium we have:

n eq = n int,ef f exp(qV /k B T ) (6)
p eq = n int,ef f exp(-qV /k B T ) (7)
Now when outside equilibrium, more specific mechanisms that must be included are:

• Mobility: we model the temperature and doping dependence of µ n and µ p via the Arora model [START_REF] Arora | Electron and hole mobilities in silicon as a function of concentration and temperature[END_REF]. We include velocity saturation (i.e. mobility degradation with the electric field) via the Canali model [START_REF] Canali | Electron and hole drift velocity measurements in silicon and their empirical relation to electric field and temperature[END_REF], and take into account mobility degradation with carrier densities (i.e. Carrier-Carrier Scattering or CCS) via the Conwell-Weisskopf model [START_REF] Choo | Theory of a forward-biased diffused-junction P-L-N rectifier -Part I: Exact numerical solutions[END_REF]- [START_REF] Fletcher | The High Current Limit for Semiconductor Junction Devices[END_REF]. Indeed, CCS has been widely recognized as one of the driving mechanisms for radiation-induced carrier transport given the very high injection levels at play [START_REF] Dodd | Device simulation of charge collection and single-event upset[END_REF]. Therefore, it needs to be included in our solver for correct behavior at high Linear Energy Transfers (LET). • Recombination: we consider Shockley-Read-Hall (SRH)

i.e. impurity-assisted recombination depending on the doping level, and Auger i.e. direct band-to-band recom-bination depending on the temperature. Auger recombination becomes especially important at high LET.

• Generation: we use the Van Overstraeten model [START_REF] Overstraeten | Measurement of the ionization rates in diffused silicon p-n junctions[END_REF] to account for impact ionization i.e. avalanche multiplication of the carriers under high electric fields. This mechanism is especially important for a proper rendition of bipolar amplification in FD-SOI, which is strongly "avalancheassisted" [START_REF] Malherbe | Alpha soft error rate of FDSOI 28 nm SRAMs: Experimental testing and simulation analysis[END_REF]. Last but not least, the radiation-induced generation of carriers is accounted for with

G rad (x, t) = N rad (x)g rad (t) (8)
where N rad (x) is the space-dependent deposited density and g rad (t) is the temporal rate at which deposition occurs. Finally, we assume electrothermal equilibrium at the contacts at all times. Using l and r subscripts for the left and right sides of the simulation domain, this leads to the following Dirichlet boundary conditions on the primary variables V , n, and p:

V (x l , t) = V eq (x l ) + V a,l (t) (9) n(x l , t) = n eq (x l ) (10)
where V a stands for the voltage applied by the circuit. Equations on the right side and for p were not written for brevity.

B. Spatial and temporal discretization

In order to yield a system of equations that can actually be solved on a computer, the set of coupled partial differential equations (PDE) above needs to be discretized on a computational grid, both in space and time. The mathematical details of the discretization procedures are given in Appendix A, here we just wish to mention the main steps taken without resorting to formulae:

• The box discretization method [START_REF] Bank | Numerical methods for semiconductor device simulation[END_REF] is used to discretize (1) (2) and ( 3), based on the general non-uniform 1D mesh of variable cross section illustrated in Fig. 1. Doing so we can account for the detailed morphology of 28 nm FD-SOI transistors, whose silicon film thickness is 7 nm under the gate, but varies near the spacers and goes to 12 nm in the epitaxy-raised source/drain implants. The electric field F is easily related to the potential V via a finite-difference formula, but the current densities J n and J p cannot be obtained from a naive discretization of ( 4) and ( 5): these equations need to be discretized using the Scharfetter-Gummel scheme [START_REF] Scharfetter | Large-signal analysis of a silicon read diode oscillator[END_REF], or else the mesh has to be extremely dense to ensure convergence. • A backward-Euler approach is used for temporal discretization, for two simple reasons: it is extremely stable (while not the most accurate), and we are anticipating on coupling with an external circuit solver: Eldo's interactive mode allows us to run a transient simulation for arbitrarily small time steps, but it is impossible to "rewind" if we want to subdivide the step further. This prohibits the implementation of a temporal discretization scheme such as the trapezoidal rule, which makes use of "two half time steps" for better accuracy. • Assembling all equations yields a non-linear system for the potentials and carrier densities, from one time step to the other. The system is solved with the Newton-Raphson method, which ultimately amounts to solving a linear system involving the Jacobian matrix of the equation residuals. LU factorization of the Jacobian is the usual procedure for 1D grids, which leads to two triangular systems solved by forward/backward sweep. • All of the above allows to propagate in time the spacedependent solution, with extrinsic boundary conditionsbut the circuit response depends itself on the device evolution. Since our circuit solver is a black box, we cannot just append the circuit equations to the device equations, and solve the entire system as one single selfconsistent entity (which is what mixed-mode TCAD often does). Therefore, additional care must be taken with the device-to-circuit communication: a basic implementation of circuit coupling is to substitute the device with a pure current source whose magnitude is calculated by the device solver, then let the circuit solver refresh the terminal voltages, and then compute a new value for the device current, and so forth. Such a coupling method was tested but exhibited mediocre stability, typically requiring time steps below 0.1 ps at all time for convergence.

Based on [START_REF] Mayaram | Coupling algorithms for mixed-level circuit and device simulation[END_REF], a much better convergence radius was then obtained by replacing the device with an equivalent current source in parallel with an equivalent conductance, as depicted in Fig. 2. Just like the current value, the conductance value is refreshed at each time step by the device solver (with negligible calculation overhead), then handed over to Eldo for circuit simulation. This allows to raise the time step to values close to the picosecond range during the impact -and much higher after -with very limited accuracy loss.

III. IMPLEMENTATION DETAILS

Our model is written in C++, in about 5,000 lines of code. During the transient run, communication with the circuit solver is implemented by redirecting the output flux of our program to the input flux of Eldo in interactive mode: using "pipe" system calls, at each time step we are able to communicate the equivalent current and conductance (device-to-circuit link), and retrieve the updated terminal potentials (circuit-to-device Fig. 2: Circuit coupling with the simulated device's equivalent conductance and current source inside the electrical solver link). Such an implementation of device-circuit connection was already presented in [START_REF] Glorieux | Random-Walk Drift-Diffusion Charge-Collection Model For Reverse-Biased Junctions Embedded in Circuits[END_REF]. Although some of the syntax is necessarily Eldo-specific, this solution can be ported to other electrical solvers as soon as they offer a commandline execution mode allowing arbitrary time increments and dynamic access to the circuit parameters. Alternatively, other works have demonstrated the use of Verilog-A models to perform a similar task [START_REF] Kauppila | A Bias-Dependent Single-Event Compact Model Implemented Into BSIM4 and a 90 nm CMOS Process Design Kit[END_REF]: these are really just different technical solutions to tackling the same challenge. Note also that before the transient simulation can be run, the solver first needs to find the equilibrium situation (zero-bias steady state, which solves (1) self-consistently with ( 6) and ( 7) starting from an initial guess based on charge neutrality), then perform a quasi-stationary ramp-up (steady state with a bias) so the device is ready to be irradiated. Now, several optimizations are needed to ensure fast execution:

• First of all, the linear system obtained for the Newton method needs to be built before it can be inverted. Calculating all of the system's protagonists is very lengthy if nothing is done to speed-up the expression evaluations. To tackle this issue, we adopt a pragmatic, hybrid strategy by caching the "clearly re-usable" arrays. For instance, the SRH recombination term is calculated with a fraction N/D on each grid point, so storing the N and D arrays allows for faster evaluation of (N D -N D )/D 2 in the Jacobian calculations later on. This requires some hard-coding of the dependence tree in the expression evaluations, to avoid refreshing values in the wrong order. Overall, we observed a near 10× speed-up when going from a careless evaluation methodology to this hybrid solution. Therefore a complete handling of expression trees -which would be required for optimal lazy reevaluation -was deemed unnecessary. • Second of all, the system needs to be inverted. For this we make use of Eigen [START_REF] Guennebaud | Eigen v3[END_REF], a very complete library for linear algebra. It allows to analyze the sparsity pattern of the Jacobian matrix to accelerate its LU factorization afterwards: the non-zero entries do vary but their indices do not, so the sparsity analysis can be performed once and for all beforehand. We observed a 40% performance Fig. 3: Comparison between our model and 3D TCAD simulations for various impacts on an NMOS inside an inverter improvement with this feature -which ironically made the system faster to invert than to build. Furthermore, Eigen enables many compiler optimizations to vectorize the calculations, yielding an additional 3× speed-up. We also implement some adaptive time step control tailored to our radiation context, with time step relaxation after the impact for instance. The end result of all these optimizations is that an impact takes about 1-3 seconds to run, depending on its toughness: this perfectly meets our speed requirements for Soft-Error Rate (SER) simulation.

IV. MODEL USAGE IN 28 NM FD-SOI

A. Validation against 3D TCAD

As a benchmark to validate our model in FD-SOI technology, we compare the results of our 1D solver with more detailed TCAD simulations. The latter are performed by striking a 3D NMOS structure connected in mixed mode to a PMOS and a load capacitance, to simulate the context of an inverter in a CMOS network. Our 1D solver simply takes as input the 1D doping profile of the struck NMOS along its channel axis. A very important note is that in our model, we do not "erase" the NMOS in the circuit simulation, unlike what was shown in Fig. 2 for a toy circuit: we cannot realistically hope for our solver to reproduce the full behavior of a MOSFET since the influence of the gate cannot be natively treated in 1D. Therefore, for the equivalent current and conductance calculations we subtract the steady-state leakage predicted by our solver, and keep the MOSFET in parallel with those circuit elements. In other words, we leave the silicon-calibrated SPICE transistor in charge of the "normaloperation" current, and our solver is only responsible for the radiation-induced current. From a device-level perspective, our simulator is accountable for the motion of excess carriers but not the "background" electrons and holes. Fig. 3 displays the current and voltage waveforms obtained for impacts at the drain-channel junction of the transistor (where the electric field is highest), for various LET. For comparison we also include results from the behavioral model in [START_REF] Malherbe | On-Orbit Upset Rate Prediction at Advanced Technology Nodes: a 28 nm FD-SOI Case Study[END_REF], where the current waveform was simply a double-exponential pulse, of integral charge given by an LET-dependent bipolar amplification modulated by a position-dependent collection efficiency. As can be seen, the time-dependent effects are well accounted for by our new model, which natively enlarges the current pulse in high injection. On the other hand with the simple model, at high LET the decoupled current source is way off since it does not vanish when the potential drops, thus causing unphysical values for the voltage. This is not an issue for Single-Event Upset (SEU) simulation which is merely interested in whether the voltage does drop or not. But Single-Event Transient (SET) simulations would be quite coarse with the simple model, since the voltage pulse duration is mispredicted. Mention should be made that the model developed in this paper is not exempt from calibration: the TCAD simulations shown here use hydrodynamic transport equations with more complex mobility modeling (with velocity saturation and interface degradation), beside the models already listed in section II-A. Thus some adjustments are necessary, mainly by playing on the carrier mobilities and lifetimes. Note however that obtaining a good match with TCAD in not an end per se: silicon measurements -e.g. comparison with SEU cross sections on sequential cells -still have to remain the main judge to close the "model calibration" loop.

B. On the importance of circuit coupling

As a mere sanity check, in Fig. 4 we plotted the waveforms obtained with our model for a 10 MeV.cm 2 /mg impact on a pn + structure with mesh dimensions representative of a bulk-Si junction as depicted in Fig. 2, with V dd = 1.2 V, R = 10 kΩ, and C = 1 fF. As can be seen, our 1D solver perfectly captures the so-called "plateau" effect [START_REF] Dasgupta | Effect of Well and Substrate Potential Modulation on Single Event Pulse Shape in Deep Submicron CMOS[END_REF], which arises in modern bulk technologies when the time constants for charge transport and circuit response start overlapping. The n + voltage is clamped at -V bi = -0.8 V, while the current lingers at the precise value of (V dd + V bi )/R = (1.2 + 0.8)/10, 000 = 0.2 mA: if the current was to exceed this value, the n + potential would drop even further and the junction would become truly forward biased, thus it would start driving a massive current in the opposite direction that would actually draw the potential back up very quickly. As a matter of fact, this phenomenon can be observed when intentionally using overly relaxed time steps during the plateau: the forward-bias current restores a dramatically high potential on the n + node, which causes spurious oscillations that sometimes stabilize, or sometimes can cause the simulation to diverge.

Voltage [V] Current [A] Time [s]

Current Voltage Fig. 4: Sanity check for circuit coupling: plateau effect for a 10 MeV.cm 2 /mg impact on a mesh representative of a bulk-Si p-n + junction connected like in Fig. 2. The plateau lasts for as long as there are charges to collect that have not recombined, and then the voltage recovers with an RC time constant.

On the other hand as can be seen in Fig. 3, in FD-SOI the plateau effect due to circuit coupling does arise, but in a less pronounced way than in bulk technology. This was expected given the reduced active silicon volumes, providing a much smaller reservoir for continued charge collection over time. More important than the current plateau however is the fact that the V ds voltage of the MOSFET is clamped at 0, much like it was clamped at V bi for a junction alone. This is of great relevance to properly capture the radiation behavior of devices connected in series. instance, consider the robust inverter shown in Fig. 5, which uses transistors of identical function in series. Stacking is a well-known method for hardening, and its efficiency has sometimes been explained by the fact that it is highly unlikely for an impact to "turn on" both off-state transistors at the same time, especially in SOI technology. More rigorously, the robustness comes from the fact that since V ds tends to clamp at 0 during an impact, for instance when we strike the NMOS closest to Z, the output cannot drop below 0.5 V for a long time and therefore, only short-lived SETs could propagate. For the voltage waveform obtained at 30 MeV.cm 2 /mg in Fig. 5, the glitch is likely to be filtered by the next logic stage. For comparison we also plotted the same simulation but without circuit coupling: the current was extracted under the same impact but with a constant 0.5 V hard bias across the NMOS, and then this frozen current waveform was injected in the circuit simulation to obtain a voltage waveform. As can be seen, in the non-coupled case the potential has no reason not to drop below 0.5 V, and therefore the glitch could very well propagate. This truly establishes the relevance of developing bias-dependent single event models if one is to try and predict the robustness of a given circuit design: without circuit coupling, the simulation simply misses the entire point of hardening by stacking in SOI. Fig. 5: In FD-SOI there is no strong current plateau, but the drain-source voltage of the struck FET does clamp at 0, highlighting the importance of circuit coupling: for a robust inverter that uses identical transistors stacked in series (a), the voltage drop can overestimated if the radiation current model does not dynamically depend on bias (30 MeV.cm 2 /mg impact).

VDD = 1 V GND = 0 V A = 0 V Z = 1 V 0.5 V (a

C. Cell-level explorations

In this last section we embed our collection model within our proprietary Monte-Carlo SER simulation platform called TIARA [START_REF] Roche | Application of the TIARA Radiation Transport Tool to Single Event Effects Simulation[END_REF]. Briefly speaking, TIARA builds a 3D model of a circuit based on its layout and technological process information. Then using the circuit's netlist, a simplified Layout Versus Schematic (LVS) is run to map the devices and nodes in the schematic to physical locations in the 3D model. The user can then specify an irradiation scenario, and ions are drawn randomly and ray-traced to compute charge deposition inside the structure. Finally, the collection model is called in conjunction with the electrical solver in order to compute the circuit's response to the impacts. Some software upgrades were needed on both sides to enable smooth integration of the new collection model within TIARA. Beside general flow unification tasks to connect the API together, we had to dynamically generate the inputs for the new collector: first, automatic generation of the 1D meshes based on the impacted transistors, and then conversion of the raytracer's output into a carrier generation term. Note that for the first task, in the case of two transistors sharing a common diffusion area we generate two simulation domains with a shared boundary condition. In Fig. 6 we plotted the results of a virtual heavy-ion irradiation on an 28 nm FD-SOI SRAM cell: experimental data was obtained at RADEF [START_REF] Virtanen | Radiation effects facility RADEF[END_REF] at room temperature, but high-temperature runs were only carried out for Single-Event Latchup (SEL) testing. Our simulations show that the saturation cross section is nearly identical at 100 • C, but at low LET it is increased by 60%. This prediction is compatible with our experimental neutron-SER measurements showing a 2× increase at 125 • C on similar bit-cells. This can be of importance for space applications, which often feature intense temperature cycles. Note that such an investigation would have been less straightforward with the simple model in [START_REF] Malherbe | On-Orbit Upset Rate Prediction at Advanced Technology Nodes: a 28 nm FD-SOI Case Study[END_REF], which is only calibrated at room temperature: it would have required a TCAD study first, then to introduce an ad-hoc temperature dependence in the collection model. In this paper we have developed a time-dependent single event model that can be interfaced with an electrical solver. The model is based on direct resolution of the semiconductor equations on a one-dimensional grid, making it computationally efficient. Furthermore this study has shown that the model is very accurate when applied to 28 nm FD-SOI devices, due to the ultra-thin silicon film in that technology. It has been integrated within our Monte-Carlo radiation simulation platform, highlighting how physics-based models can natively enable advanced sensitivity analyses at cell level. Note that the model could be applied seamlessly to more advanced nodes in FD-SOI -e.g. 22 or 12 nm -and, with some adaptation, to FinFET technology as well: using simplified 3D diffusion in the bulk beneath the fin and looking at carriers collected at the base of the fin, would allow to derive an effective generation term inside the fin and then, use 1D transport as shown in this article. Note that this implicitly assumes that collection inside the fin does not strongly affect transport in the bulk beneath, which seems reasonable since the electric field is mostly confined within the channel. From a physicist's perspective, compact analytical models are often preferred for the immediate insight they can offer, while brute-force numerical approaches are sometimes regarded as merely "number-crunching" routines making it hard to build one's physical intuition. However in this case we find that, owing to the solver's great speed, having the ability to get an instantaneous answer to any thought experiment -in 1D -is also extremely valuable and helps us gain a deeper understanding of the underlying mechanisms.

APPENDIX

Our goal here is to shed more light onto the mathematical details of the discretization procedure that was described succinctly in Section II-B. Suppose we want to compute the solution to our coupled PDE set on discrete grid points {x i } i=0..I-1 and discrete instants {t j } j=0..J-1 . We will note f (x i , t j ) = f i,j for short. We also scale our units to work on dimensionless quantities: all potentials are divided by the thermal voltage V T = k B T /q, all densities (carriers, doping, radiation...) are divided by the intrinsic density n int , and to close the loop, all times are divided by an arbitrary, extrinsic, reference time. First, the box method [START_REF] Bank | Numerical methods for semiconductor device simulation[END_REF] is used for spatial discretization; for all equations of the form ∇ • F + S = 0 i.e. (1) (2) and (3), the idea is to use Ostrogradky's theorem on a local control volume (the "box" depicted in Fig. 1), yielding:

σ i+ F i+ -σ i-F i-+ v i S i = 0
In Poisson's equation for instance, this replaces the usual finite-difference formula for ∆V with a reweighted Laplacian operator based on the local mesh cross section, natively leading to steeper potential drops in areas where the silicon is thinner. More precisely, if we denote dx 2 i-= v i (x i -x i-1 )/σ i- and dx 2 i+ = v i (x i+1 -x i )/σ i+ , then the Laplacian stencil is expressed as:

∆ i = 1 dx 2 i- , - 1 dx 2 i- + 1 dx 2 i+ , 1 dx 2 i+
Now the current equations ( 4) and ( 5) are more specific, and need to be discretized via the Scharfetter-Gummel scheme [START_REF] Scharfetter | Large-signal analysis of a silicon read diode oscillator[END_REF]. Briefly summarized, since the electron and hole densities have exponential dependence on the potential, when using naive finite-difference formulae the equations are made locally inconsistent unless the mesh is extremely dense. The Scharfetter-Gummel method prevents this issue by anticipating the values of n and p in-between grid points, yielding e.g. for the leftward current densities:

J n,i-= µ n,i-[n i B(V n,i -V n,i-1 ) -n i-1 B(V n,i-1 -V n,i )] J p,i-= µ p,i-[p i-1 B(V p,i -V n,i-1 ) -p i B(V p,i-1 -V p,i )]
where the Bernoulli function:

B(x) = x e x -1
performs an optimal transition from a usual finite difference formula (which works well for diffusion-driven transport) to an upwind difference formula (which works well for driftdominated problems) depending on carrier velocity. Note that B(x) needs to be evaluated carefully to avoid overflows at large values of x and imprecision due to the indeterminate form when x → 0 (in that case we use a Taylor series expansion evaluated to machine precision). As for temporal discretization, with the backward-Euler method all expressions of the form ∂g/∂t+h = 0 are replaced by (g j -g j-1 )/(t j -t j-1 ) + h j = 0. Then when going from j -1 to j, assembling all discretized equations (with (4) and (5) substituted into (2) and ( 3)) yields a non-linear system of 3I equations for 3I unknowns -which we shall refrain from writing here in its entirety for the sake of compactness; if we denote our state vector X, then it is of the form: 0 3I = f (V 0,j ..V I-1,j , n 0,j ..n I-1,j , p 0,j ..p I-1,j ) = f (X j )

The root to the equation residuals f is computed via the Newton-Raphson method, which performs successive refinements to the solution by linearizing the system. To go from "outer" iteration k to k + 1 we compute an update candidate dX k+1 j as the solution to the linearized system:

0 3I = J f • dX k+1 j + f (X k j )
where J f stands for the Jacobian matrix of f . Then to actually refresh the solution we compute X k+1 j = X k j + λ k+1 dX k+1 j where λ is a damping parameter taken "as close as possible to 1", based on whether or not the update actually does improve the residuals -typically with geometrical damping for the "inner" iterations. Overall convergence is checked by ensuring the norm of the right-hand side falls below a given value.

Optionally, we could also use a criterion the norm of the solution update, but since the Newton method converges quadratically (doubling the number of exact digits at each step) when initialized close enough to the solution, this could cause to systematically perform one useless iteration just to let the solution stabilize when it is already good enough. A last remaining mathematical aspect is the calculation of the equivalent conductance for our simulated device. It is defined as G eq = ∂I/∂U , where the applied voltage difference is U = V l -V r -(V 0,eq -V I-1,eq ) and current I is obtained via the mesh average of the flux of the current densities in order to minimize numerical error:

I = 1 I -1 I-2 i=0
(J n,i+ + J p,i+ ) • σ i+ (the displacement current I d can safely be ignored: even when the voltage drops by 1 V in 1 ps, for our L ≈ 100 nm, σ ≈ 10 × 100 nm 2 meshes it only reaches I d ≈ εσ/L • ∂U/∂t ≈ 1 µA). Intuitively, the derivative for G eq should not have to be computed by evaluating the solution "at U + dU " -which could double the overall computation time -since the Jacobian matrix J f already contains all the information about how the equation residuals vary around the solution. This is precisely the idea in [START_REF] Mayaram | Coupling algorithms for mixed-level circuit and device simulation[END_REF]: we can write

G eq = ∂I ∂X
• ∂X ∂U The first term is evaluated from symbolic differentiation, and to obtain the second term we make use of the fact that f (X(U ), U ) = 0 =⇒ J f • ∂X ∂U = -∂f ∂U ∂f /∂U is a nearly-empty vector since U only appears at the boundary nodes, and because the Jacobian is already LUfactorized, the above system can be solved for ∂X/∂U at almost no computational overcost.

Fig. 1 :

 1 Fig. 1: Notations for 1D mesh of variable thickness along x

  -10 1.01E-09 1.03E-09 1.05E-09 1.07E-09 1.09E-09

Fig. 6 :

 6 Fig. 6: SRAM cross section simulated at elevated temperature


Victor Malherbe

Gilles Gasiot

Dimitri Soussan

Jean-Luc Autran

Philippe Roche

On-Orbit Upset Rate Prediction at Advanced Technology Nodes: a 28 nm FD-SOI Case Study

We address accurate computation of on-orbit upset rates in advanced technologies, with a focus on FD-SOI at the 28 nm node. Heavy-ion measurements performed on FD-SOI SRAM bit-cells give experimental evidence of the technology's intrinsic robustness in space environments; this extreme reduction of sensitive volume dimensions deeply affects the assumptions pertaining to the radiation response models used to predict upset rates. The generic "Integral Rectangular ParallelePiped" (IRPP) model, although requiring careful setting of its parameters, provides us with first-order estimates of the error rate. We then present a custom FD-SOI response model within our Monte-Carlo simulation chain, enabling comparison with IRPP and further analyses.

I. INTRODUCTION

Accurate prediction of the rate of Single-Event Upsets (SEU) is of utmost importance for space missions featuring a harsh radiation environment: while under-prediction of the system's SEU rate may result in mission failure, overestimating it can lead to an overly conservative design coming at a high cost, especially given the timescales involved. Because direct measurements of SEU's on orbit are quite rare, the most practical way to enable projections of a circuit's error rate in space is to use a combination of accelerated test data along with space environment models describing the particle spectra encountered on orbit. In many cases, the irradiation experiments can only be performed at normal incidence or just a few angles due to limited beam time; therefore, an extra modeling effort is required to extrapolate the response of the technologies under an omnidirectional space environment. To that end, several SEU rate prediction methodologies have been developed over the years, for instance the Integral Rectangular ParallelePiped (IRPP) method originally developed in [START_REF] Petersen | Rate prediction for single event effects-a critique[END_REF] and now implemented in popular tools such as CREME [2] and SPENVIS [3]. However, emerging technologies such as Fully-Depleted SOI (FD-SOI) or FinFET feature strongly-confined geometries that are challenging the existing SEU models [START_REF] Roche | Advanced CMOS Bulk, FinFET and UTBB SOI Technologies[END_REF]. In this study we thus wish to discuss what modeling approaches are relevant in advanced depleted technologies. This paper is organized as follows: in sections II and III, experimental heavy-ion cross sections are presented for a 28 nm FD-SOI test vehicle embedding SRAM arrays, and compared with measurements on bulk CMOS counterparts. The next two sections IV and V discuss how these measurements are used as a basis for extrapolation to real space environments: first-order predictions of the SEU rate are obtained via the IRPP method, and we examine the validity of some crucial aspects in the model. Then, by means of Monte-Carlo simulations tightly calibrated with knowledge of the technology's fabrication process and radiation response, we derive a best estimate for the error rate in our study case. This provides us with generic recommendations regarding the usage of CREME for advanced technologies featuring narrow sensitive volumes such as FD-SOI or FinFET.

II. TEST CHIP AND EXPERIMENTAL SETUP

The irradiated test vehicle was designed and manufactured in STMicroelectronics' 28 nm FD-SOI technology. GDS and real views of the layout are given in Fig. 1: the test chip embeds several flip-flop registers and SRAM cuts, as well as microprocessors with various implementations (reference, radiation hardened...). The heavy-ion testings were performed at RADEF, Finland [START_REF] Virtanen | Radiation effects facility RADEF[END_REF], in compliance with ESA/SCC basic specification No 25100. Five ion cocktails were used, namely nitrogen, neon, iron, krypton and xenon ions giving rise to Linear Energy Transfers (LET) of 1.83, 10.2, 18.5, 30.4 and 60.0 MeV.cm 2 /mg. Three parts were irradiated with fluences of about 10 7 ions/cm 2 on each run and, except for the Single Event Latchup runs, room temperature conditions were used and the memory arrays' supply voltage was V dd,mem = 1.0 V.

III. HEAVY-ION TEST RESULTS

The overall outcome of this test campaign is a very high intrinsic resilience exhibited by the FD-SOI technology under heavy ions. Previous publications had already reported the outstanding SEU hardness of 28 nm FD-SOI under terrestrial particles such as alphas and neutrons [6]- [START_REF] Gasiot | SER/SEL performances of SRAMs in UTBB FDSOI28 and comparisons with PDSOI and BULK counterparts[END_REF]; this work confirms the technology's hardness in spatial environments as well, with heavy-ion cross sections routinely achieving 2 decades of reduction when compared to equivalent-node bulk CMOS technologies. As far as SEU's are concerned, this extreme resilience of FD-SOI is primarily due to two factors: first of all the very small volume for ion-induced charge deposition, thanks to the 7 nm-thin active silicon film being enclosed by the Buried Oxide (BOX). Second of all, parasitic bipolar amplification is quite limited in this technology, as was investigated with TCAD simulations in [START_REF] Malherbe | Alpha soft error rate of FDSOI 28 nm SRAMs: Experimental testing and simulation analysis[END_REF], leading to collected charges as small as 0.1 to only a few fC. To help narrow down our discussion, the rest of this study will focus on just one of the SRAM cells tested. The bit-cell follows a classical design with six transistors (Fig. 2), out of which three are SEU-sensitive depending on whether a 0 or a 1 is stored. For instance, in retention mode the bit-and wordlines are not loaded (BLT = BLF = W LL = W LR = 0) and if the cell is at 1 (BLT I = BLF I = 1), then only PD1 (pull-down NMOS), PG1 (pass-gate NMOS) and PU2 (pullup PMOS) will host an electric field sufficient to separate the radiation-induced electron-hole pairs. The experimental heavy-ion cross sections collected in the beam tests for our FD-SOI bit-cell are shown in Fig. 3 along with a Weibull fit. As an element of comparison, the cross section for a 28 nm "bulk-equivalent" bit-cell was also plotted (identical area and similar design, give or take some Design Rule Check alterations). As can be seen, the FD-SOI cell is about two decades less sensitive than the bulk one: the bulk bitcell exhibits an asymptotic cross-section of about 20 times the cell area, meaning that on average each impact triggers a 20-cell upset, because the carriers are free to diffuse on very long ranges. On the other hand, the FD-SOI limiting cross section σ sat is only about 15% of the cell area thanks to the full dielectric isolation provided by the BOX. In bulk technology, the primary sensitive area is that of the drains of blocked transistors (with a preponderant vertical field) and it only requires an LET of about 1.5 MeV.cm 2 /mg to be "activated". However in FD-SOI the situation is very different: the sensitive area is rather located by the off-transistor channels featuring a strong lateral field -and it takes an LET of about 15 MeV.cm 2 /mg to "activate" it entirely. Such considerations will be used to construct the sensitive volumes for charge collection in the next sections where error rates are derived. The aim of this section is to describe the flow followed to compute on-orbit SEU rates in FD-SOI using CREME tools provided online [2]- [START_REF] Tylka | CREME96: A Revision of the Cosmic Ray Effects on Micro-Electronics Code[END_REF]. We first simulate ion environments encountered under various orbit and solar conditions, and then make use of the IRPP method to convolve the environment with the experimental test data. We also discuss the influence of crucial parameters in the IRPP setup.

A. Radiation environment simulation

Several radiation environments were simulated in CREME, ranging from Low Earth Orbit (LEO) to Geosynchronous Earth Orbit (GEO). Trapped protons are accounted for with the AP8 model, while Galactic Cosmic Rays (GCR) are calculated with CREME96 GCR model. We cover multiple solar conditions, namely minimum and maximum "solar quiet" configurations (respectively maximizing and minimizing GCR fluxes) and "worst week" solar flare. All ion species up to Z=92 (uranium) are considered for flux calculations, and transported through 100 mils of aluminium. A sample of resulting LET spectra is plotted in Fig. 4, highlighting differences not only in total flux magnitude, but also in LET content: the GEO environment under quiet Sun contains heavier species, while the LEO or flare environments are dominated by light particles. Note that similar fluxes were obtained with SPENVIS and are not reported: the scope of this study is not to discuss the accuracy of radiation environment models, but rather on the SEU models taking those fluxes as an input.

B. IRPP method and discussion

The IRPP method [START_REF] Petersen | Rate prediction for single event effects-a critique[END_REF], due to its genericity and large case coverage, has been widely adopted in the radiation effects community to compute SEU rates in space from groundbased test data. First of all, given a fixed critical charge Q crit , the (non integral) RPP method works by considering a sensitive volume of x, y, z dimensions chosen to match both the experimental saturation cross section at normal incidence: x • y = σ sat and threshold LET:

LET th • z = Q crit .
Then the angular response is extrapolated by assuming that the sensitive volume has uniform 100% collection efficacy: for an ion of given LET and direction, this means that the collected charge will equal the charge deposited within the sensitive volume, yielding Q coll = Q dep = LET • l where l is the chord length traveled by the ion within the sensitive volume. This collected charge gives rise to an upset when it is superior to Q crit . Thus, in the RPP method, the predicted cross section at normal incidence is necessarily a step function. The more general IRPP method takes the above ideas further, by interpreting the σ(LET )/σ sat relative cross-section curve as a cumulated distribution function for the threshold LET's of RPP's of identical dimensions: while a step-function cross section curve corresponds to a single threshold LET, a gradual Weibull curve starting with an onset LET represents a wider distribution of threshold LET's from the onset LET value onwards. Note however that the physical meaning of a statistical superposition of identical sensitive volumes with different threshold LET's is quite questionable: at the time in 1992, it was believed that each individual bit-cell would exhibit a sharp cross-section curve, and that the smooth curve measured was due to a spread across the memory array due to process variability. We now Fig. 5: GEO upset rates obtained with CREME via the IRPP model for FD-SOI and bulk technologies (a) for various RPP dimensions and (b) for various solar activities know that even one cell alone intrinsically exhibits a gradual activation, which is the motivation for using "nested" sensitive volumes in certain Monte-Carlo codes such as MRED [START_REF] Reed | Physical Processes and Applications of the Monte Carlo Radiative Energy Deposition (MRED) Code[END_REF]. Given an experimental cross section curve, the x and y dimensions of the IRPP can reasonably be estimated, starting from √ σ sat (the default setting in CREME) and then playing with various aspect ratios. However, the choice of the z dimension can be quite blurry; physically, it represents the depth over which quick collection by electric drift occurs, i.e. the depletion depth. In bulk CMOS technologies, the usual recommendation in CREME is to set z around 1 µm. This is unrealistic for our 28 nm FD-SOI technology here, whose active silicon thickness is 7 nm (charges deposited in the BOX or below are not collected). SPICE simulations were performed on the studied bit-cell and revealed critical charges Q crit of about 0.4 fC (resp. 0.8 fC) for NMOS (resp. PMOS) impacts.

Together with LET th = 1.5 MeV.cm 2 /mg ⇐⇒ 15.45 fC/µm, this suggests an effective collection depth of a few dozen nanometers, which is already a lot more in line with the physical sensitive volume. Fig. 5 (a) depicts the z dependence of the SEU rate predicted in GEO for the FD-SOI and bulk bit-cells with the IRPP method, using the Weibull fits from Note that the increase in SEU rate with decreasing z, although counter-intuitive at first ("how can a smaller volume be more sensitive?"), is easily understood either in geometrical terms (few chords longer than z for a "needle", and the opposite for a "slab") or in electrical terms (low Q crit for a small z). The IRPP x:y aspect ratio was also investigated using information from the cell layout: since the NMOS critical charge is much lower than for the PMOS and since they have much larger area, the sensitive volume x dimension can reasonably be approximated by the NMOS width. However, as shown in Fig. 5 (a), this has very little effect on the SEU rate prediction, giving no more than 5% correction from the "square IRPP" predictions -and the difference becomes negligible as z gets much smaller than x and y, i.e. in the "slab" configuration where chord lengths are driven by z.

Using z = 10 nm for FD-SOI and z = 1 µm for bulk, in Fig. 5 (b) we compare both technologies' upset rates for various conditions. As can be seen, it takes a "worst week" solar flare for FD-SOI to suffer more upsets than bulk in quiet solar conditions. In Fig. 6 we also study the influence of orbit parameters: for the LEO's, the sharp difference between 800 km and 2000 km reflects the presence of the first Van Allen belt, with higher trapped proton populations above approximately 1000 km. The more modest increase in upset rate from 15000 km to 30000 km is not due to trapped particles (since the second Van Allen belt is populated with electrons that are mostly stopped by 100 mils of shielding), but rather by the depletion of the Earth's magnetic shield towards higher altitudes. Finally, at altitudes as low as 800 km the orbit inclination can have a strong impact, especially in the South Atlantic Anomaly.

Note that throughout this study we only consider direct ionization by the ion species on orbit, thereby neglecting upsets from secondary nuclear products: we do account for the trapped protons' contribution to the received LET spectrum, but overlook the upsets triggered by ionizing products of protons against silicon. At LEO altitudes (2000 km and below) such a choice might not be justified, but it allows comparison on a common ground for the various models under discussion.

V. BEST ESTIMATE OF THE SEU RATE IN TIARA

In this section we go through the derivation of a more accurate SEU rate by means of a dedicated response model for FD-SOI. We compare the results obtained with those from CREME and give insights into their physical meaning.

A. Overview of TIARA flow

TIARA stands for ToolsuIte for rAdiation Reliability Assessment, and is a proprietary tool developed at STMicroelectronics to study Single Event Effects (SEE) at circuit level [START_REF] Roche | Application of the TIARA Radiation Transport Tool to Single Event Effects Simulation[END_REF]. TIARA's flow for prediction of SEE's can be summarized followingly:

• Structure generation from circuit netlist and layout: using the x and y coordinates from the GDS together with process information ("the z coordinates" such as implantation depths, Back-End-Of-Line or BEOL dielectric thicknesses), a 3D structure is generated for the circuit.

B. FD-SOI radiation response model and calibration at normal incidence

The main specificities of our 28 nm FD-SOI technology under radiations have been investigated with TCAD simulations and were reported in [START_REF] Malherbe | Alpha soft error rate of FDSOI 28 nm SRAMs: Experimental testing and simulation analysis[END_REF]. The primary parameter to account for, beside the sensitive volume small dimensions, is the bipolar amplification defined as the ratio of the collected charge to the deposited charge. In an NMOS for instance, radiation-induced holes in the body region tend to lower the potential barrier and trigger the parasitic lateral bipolar transistor, giving rise to current injection by the source. This phenomenon is intrinsic to all SOI technologies where majority carriers can linger in the body region. In 28 nm FD-SOI, the β factor was shown to be no more than 4 for LET's above 1 MeV.cm 2 /mg. Expanding on the work done in [START_REF] Malherbe | Alpha soft error rate of FDSOI 28 nm SRAMs: Experimental testing and simulation analysis[END_REF] where a constant bipolar amplification was accounted for in TIARA, we devise our response model for FD-SOI with a spatially-varying collection efficiency η to modulate the LET-dependent bipolar amplification along the gate axis of each sensitive transistor. To capture the non-uniformities of this charge collection "yield", the ion track inside the active silicon is discretized at points y i (with nanometric resolution typically). Then, the total collected charge is computed as:

Q coll = i Q dep (y i ) • β(LET ) • η(y i ) (1)
Using this method, for any impact we are able to calculate the collected charge Q coll . It is then related to a doubleexponential current waveform of rise and fall times τ r and τ f by:

i(t) = Q coll τ f -τ r • exp(-t/τ f ) -exp(-t/τ r ) (2)
and this is what ultimately feeds the SPICE solver. As a summary, the "free" parameters in our response model are the LET-dependent bipolar amplification, the positiondependent yield, and the current waveform rise and fall times.

To fix all of these, we use TCAD simulations to provide us with orders of magnitude and relative trends, while still relying on calibration against heavy-ion data at normal incidence to set the precise figures:

• β(LET ) relative trends for NMOS and PMOS were obtained in TCAD, and precise numbers are obtained to match the observed threshold LET (Fig. 8 (a)). • η(x) is a non-uniform weight between 0 and 1 along the gate axis (Fig. 8 (b)). In accordance with works such as [START_REF] Castellani-Coulie | Simulation analysis of the bipolar amplification in fully-depleted SOI technologies under heavy-ion irradiations[END_REF], it is set at 1 at the drain-channel junction and goes to 0 further away. Note that although the drain is strongly doped (while the channel is intrinsic), charge collection does not vanish immediately inside it. Defining the spatial dependence of this yield truly requires nanometric resolution around the junction, for the experimental heavyion cross section to be matched. (Fig. 8 (c)). Note that for other gate lengths than L min , the function is simply stretched by the L g /L min ratio. This is confirmed by our TCAD simulations, showing only small variations in collection for other transistor lengths. • τ r and τ f the rise and fall times of the double-exponential waveforms are both below 10 ps, representative of the current shapes computed in TCAD; since the critical charge, roughly speaking, scales as C.V +I.τ , modulating such small values of rise and fall times with the LET or the position would make little difference to the end result. • Finally, further simulations suggest that the supply voltage V dd only has a weak influence on all of the above parameters (e.g. 10% decrease in Q coll for -50% V dd ): while in bulk technology, the depletion depth strongly expands with V dd , resulting in deeper collection, in depleted technologies like FD-SOI the collected charge is primarily dictated by the -rigid -sensitive volume. Using those settings for our parameters, we reproduce the heavy-ion measurements with an excellent agreement. We then perform some sanity checks to ensure that the angular response is properly captured. Namely, simulations of alpha Soft Error Rate (SER) were performed and confronted to experimental test results across a variety of bit-cells and supply voltages. The α-SER was predicted with 1.5× accuracy under

C. TIARA SEU rate calculations in FD-SOI

Once the calibration phase is performed, the response model in TIARA can be used to predict upset rates in arbitrary environments. Provided an isotropic integral LET spectrum given by CREME (φ(LET ) in particles per unit time, area and steradian, with total flux φ(0)), the most direct way of calculating the upset rate is to consider the normalized flux φ(LET )/φ(0) to be a probability distribution for drawing LET's randomly in the Monte-Carlo simulator. This reproduces the natural distribution in the environment (alternatively, importance sampling could be implemented but is not discussed here). Under those conditions, the upset rate R in events per unit time is simply given by:

R = 4π • φ(0) • A irrad 4 • N SEU N impacts (3)
where N SEU is the recorded upset count, N impacts is the total number of particle strikes, and A irrad is the area of the irradiated surface which must define a volume enclosing all sensitive volumes (in TIARA we use either a bounding sphere or a bounding box). The 4π factor accounts for the total solid angle, and the A irrad /4 term arises as the average projected area of the irradiation surface when viewed from all angles (as long as it is convex). On a computational note, let us mention that some acceleration is needed for the simulations to run within a reasonable time: for sufficient statistics to be gathered, a few million impacts have to be simulated for each environment, which even with parallelization requires too much computing power. To avoid running all impacts entirely, we discard the calculations if the LET is below a certain cut-off, calculated to be the absolute minimum for an SEU to occur (by considering the longest chord in our sensitive volumes, the maximum bipolar amplification, the minimum critical charge...). For all environments, a huge portion of the LET spectrum stands below this cut-off, which provides a speed-up factor of about 100, resulting in simulation times of about a few hours for one environment. Direct comparisons between the IRPP upset rates computed in section IV and those calculated with TIARA are given in Fig. 9 (showing error bars from TIARA was deemed unnecessary here since all simulations gathered more than 400 SEU's, with 95% confidence intervals of ±10% = 2/ √ 400). The IRPP model provides satisfactory estimates of the error rates when TIARA is taken as a reference: overall it correlates with about 2× precision to our custom response model, with a mismatch depending on the environment LET content, as was to be expected from the LET-dependent modeling of bipolar amplification. Although both methods are in perfect agreement at normal incidence by construct, the notable difference is that with the IRPP model the gradual turn-on of the sensitive volume is forced via a statistical superposition of threshold LET's, when in our specific response model it is explained from the physical insight gained with TCAD. The legend of Fig. 9 aims at depicting this ground difference, showing nested sensitive volumes for TIARA and stacked RPP's for IRPP (even though both models simulate a continuum of those volumes in fact). Nevertheless, for practical purposes this study shows that first-order estimations of the upset rates can be obtained with the IRPP algorithm when taking the silicon film thickness as IRPP z depth. This is mainly due to the fact that parasitic bipolar amplification remains modest at all times in FD-SOI; in Partially-Depleted SOI technologies often featuring β larger than 10 [START_REF] Ferlet-Cavrois | Characterization of the parasitic bipolar amplification in SOI technologies submitted to transient irradiation[END_REF], this conclusion would probably not hold as the effective collection depth could be much larger than the physical silicon thickness above the BOX. Taking advantage of the predictive power of the circuit solver coupled to our radiation response model, we can perform further analyses with TIARA: in Fig. 10 we show the effect of voltage scaling on the technology radiation response for different orbits. The simulations show that the upset rate can increase twofold when lowering the voltage by 20%. The reason for this was mentioned earlier: while the critical charge scales linearly with voltage, the collected charge does not vary much with voltage in depleted technologies, and thus the balance between them can shift significantly. Note that the effect is symmetrical: strong reductions in SEU rate can also be achieved by raising the supply voltage, providing interesting opportunities for hardened designs. Fig. 10 also points out that the voltage scaling effects are more pronounced in low-LET environments where many ions are close to threshold LET. As an indication, the simulated upset rate for FD-SOI in GEO at solar minimum is below 5 SEU/Gb/day at all V dd , about 50 times lower than bulk simulated with IRPP.

VI. CONCLUSION

In this study we have discussed methodologies for computing on-orbit Single-Event Upset (SEU) rates for a 28 nm FD-SOI technology. Heavy-ion cross sections of a representative SRAM cell were reported to be two decades lower than an equivalent cell in planar bulk technology. Using simulations of several on-orbit radiation environments, these test data were used to compute upset rates with the Integral Rectangular ParallelePiped (IRPP) model available CREME, once again highlighting fundamental differences with bulk technology: for reduced sensitive volumes in technologies such as FD-SOI or FinFET, confinement deeply impacts the modeling choices and assumptions to be made regarding the cell's radiation response at an angle, given measurements at normal incidence: on our FD-SOI example the upset rate was shown to increase by 30× when varying the sensitive thickness in the IRPP model from the CREME default setting for bulk to the physical film thickness. For practical usage of the IRPP model in CREME, this study has shown that in 28 nm FD-SOI using the silicon film thickness as sensitive depth is a good choice, since bipolar amplification is limited in the technology. More accurate upset rates were then computed within the TIARA Monte-Carlo tool chain via a new response model tailored to FD-SOI, after TCAD investigations and calibration against heavy-ion measurements. The SEU rates calculated with IRPP were shown to be in good (2×) correlation with the reference TIARA simulations. Additionally, the effects of voltage scaling were investigated with TIARA, revealing 2× variations in upset rate when moving 20% away from nominal voltage, as a very consequence of the technology being fully depleted. For geosynchronous earth orbits under quiet Sun conditions, the SEU rate in 28 nm FD-SOI was found to be of only a few SEU per Gb per day (about 50× better than in bulk), demonstrating the technology's suitability for space applications.

Fig. 1 :

 1 Fig. 1: Chip layout (a) -photograph after package etching (b). The SRAM cuts occupy the left half of the design.

Fig. 2 :

 2 Fig. 2: The studied 6-transistor SRAM bit-cell: (a) its netlist and (b) its -idealized -layout

Fig. 3 :

 3 Fig. 3: Experimental heavy-ion cross sections of the bit-cells, exhibiting FD-SOI's intrinsic robustness compared to bulk

Fig. 4 :

 4 Fig. 4: Integral LET spectra simulated with CREME in various environments: GEO under quiet Sun features high-LET species, while LEO or active Sun environments are dominated by low-LET protons and alphas

Fig. 6 :

 6 Fig. 6: Upset rates obtained with CREME via the IRPP model for FD-SOI (a) for various altitudes at 0 inclination and (b) for various inclinations at 800 km

Fig. 7 Fig. 7 :

 77 Fig. 7: TIARA 3D structure for FD-SOI bit-cell simulation: (a) Global view up to the power rails (not displayed for clarity sake: all dielectrics except the BOX, substrate underneath shallow trench isolation) (b) Close-up on alpha particle tracks triggering SEU's within sensitive NMOS channels

Fig. 8 :

 8 Fig. 8: Breakdown of TIARA's dedicated response model for FD-SOI: (a) LET-dependent bipolar amplification -(b) Nonuniform collection efficiency -(c) Simulated heavy-ion cross section after calibration phase

Fig. 9 :

 9 Fig. 9: Comparison of TIARA and IRPP upset rate predictions under different environments: the generic IRPP approach can correlate to our custom response model with 2× accuracy.

Fig. 10 :

 10 Fig. 10: Investigation of voltage scaling effects with TIARA: the upset rate can increase by 2× upon -20% reduction in supply voltage.