Multi-timescale microscopic theory for radiation degradation of electronic and optoelectronic devices

A multi-timescale hybrid model is proposed to study microscopically the degraded performance of electronic devices, covering three individual stages of radiation effects studies, including ultrafast displacement cascade, intermediate defect stabilization and cluster formation, as well as slow defect reaction and migration. Realistic interatomic potentials are employed in molecular-dynamics calculations for the first two stages up to 100\,ns as well as for the system composed of layers with thickness of hundreds times of lattice constant. These quasi-steady-state results for individual layers are input into a rate-diffusion theory as initial conditions to calculate the steady-state distribution of point defects in a mesoscopic-scale layered-structure system, including planar biased dislocation loops and spherical neutral voids, on a much longer time scale. Assisted by the density-functional theory for specifying electronic properties of point defects, the resulting spatial distributions of these defects and clusters are taken into account in studying the degradation of electronic and optoelectronic devices, e.g., carrier momentum-relaxation time, defect-mediated non-radiative recombination, defected-assisted tunneling of electrons and defect or charged-defect Raman scattering as well. Such theoretical studies are expected to be crucial in fully understanding the physical mechanism for identifying defect species, performance degradations in field-effect transistors, photodetectors, light-emitting diodes and solar cells, and in the development of effective mitigation methods during their microscopic structure design stages.


Introduction
Point defects (vacancies and interstitial atoms) are produced by the displacements of atoms from their lattice sites, Sigmund, 2006) where the atom displacements are mainly induced by a Primary Knockout Atom (PKA) on a time scale shorter than 50 ps. This initial phase is followed subsequently by a defect reaction (clustering or dissolution of clusters), (Devanathan et al., 2001) and further by the thermally activated migration (Posselt et al., 2005) of the point defects and defect clusters over a time scale longer than 100 ns. The combination of all these processes, resulting in a significant concentration of surviving defects in the crystal, is physically termed particle irradiation displacement damage (in addition to the well know γ-ray electron ionization cascade damage). Such radiation displacement damage effects depend not only on the energy-dependent flux of the incident particles (protons, neutrons, ions, etc.) but also on the differential energy transfer cross sections (probabilities) for collision between atoms, interatomic Coulomb interactions and kinetic-energy loss to electrons inside an atom. Irradiation temperature also significantly affects the motion of defects, their stability as clusters and the formation of Frenkel pairs (Gao and Weber, 2003a).
On the other hand, electron devices are usually classified either as electronic ones, where electrons respond to an applied voltage as a current flow, or optoelectronic ones, in which electrons perform interband/intraband optical transitions in the presence of an incident light signal. (Gumbs and Huang, 2013) For an electronic device (e.g., field-effect transistors in an integrated circuit), the momentum-relaxation time of electrons, due to scattering by randomly-distributed defects, plays a crucial role in determining the electron mobility, (Strour et al., 2003) while the photo-excited electron lifetime, due to nonradiative recombination with defects, is proven to be a key factor affecting the sensitivity or the performance of optoelectronic devices (e.g., photo-detectors and lightemitting diodes) (Weatherford and Anderson, 2003).
In a perfect crystal, the continuous free-electron states are quantized into many Bloch bands separated by energy gaps and these Bloch electrons move freely inside the crystal with an effective mass different from that of the free electrons (Callaway, 1991). In the presence of defects, however, the field-driven current flow of Bloch electrons in a perfect crystal will be scattered locally by these defects, leading to a reduced electron mobility. In addition, photo-excited Bloch electrons could acquire a shortened lifetime, giving rise to a degraded quantum efficiency due to enhanced non-radiative recombination with defects. The dangling bonds attached to the point defects may capture extra electrons to form charged defects. In this case, the positively-charged holes in the system will be trapped to produce a strong spacecharge field, while the negatively-charged electrons may generate the so-called 1/f-current noise in their chaotic motion due to the presence of many potential minima and maxima from randomly-distributed charged defects.
Point defects in crystals, as shown in Fig. 1, can be generated by particle irradiation in both bulk and nanocrystals composed of many grains with different sizes . One of the effective calculation methods for studying the non-thermal spatial-temporal distributions of radiation-induced point defects is the Molecular-Dynamics (MD) model based on a stepped time-evolution approach (also termed the collisional and thermal spike stages), which involves the total force by summing over the interatomic potentials from all the atoms in a finite system (Gao et al., 2009). The lattice vibration at finite temperatures can be taken into account by an initial thermal-equilibrium state for atoms (intrinsic vacancies and interstitials) in the system plus an initial velocity for one of the atoms in a specific direction. The system size increases quadratically with the initial kinetic energy of particles and the time scale runs up to several hundred picoseconds (called the quenching stage). Therefore, the defect reaction process by thermal migration cannot be included in this MD model due to its much longer time scale, although the other processes, such as displacement of lattice atoms, energy dissipation, spontaneous recombination and clustering, can be fully taken into account. If the system time evolution goes beyond 100 ps, the kinetic lattice Monte-Carlo method can also be used (Rong et al., 2007). However, if the time scale exceeds several hundreds of nanoseconds (also called the annealing stage), the rate theory (Maksimov and Ryazanov, 1980;Golubov et al., 2012) has to be called in for studying the steady-state properties of the surviving defects (up to hours or days or even months).
We know that the MD model with a realistic interatomic potential has been developed for studying the nonthermal spatial-temporal distributions of radiation-induced point defects in noble transition metals and alloys and the density-functional theory has been widely used for calculating electronic properties of defects with pre-assumed specific defect configurations. On the other hand, a quantummechanical model has been well established for investigating defect effects on semiconductor electronic devices in the presence of spatially-uniform and randomly distributed point defects. However, to the best of our knowledge, no first-principle model and theory has been proposed so far to study microscopically the degraded performance of electronic devices induced by particle irradiation displacement damage.
Therefore, the theory presented in this study is expected to be very important in understanding the full mechanism for characterizing defects, performance degradations in transistors, photo-detectors, lightemitting diodes and solar cells, as well as in developing effective mitigation in early design stages. Equipped with our current multi-timescale microscopic theory, the experimental characterization of post-irradiated test devices will result in useful information on the device architecture's susceptibility to space radiation effects (we already know details on particle irradiation sources). On the other hand, the model should allow prediction of devices degradation on the basis of space weather forecast for the particular orbit.
Some of the equations presented below will be well-known to researchers in the materials science field, however, researchers in the device physics field may not be aware of them. With this paper, we hope to bridge the gap between researchers studying radiation-induced damage in materials and researchers studying radiation-induced performance degradation in devices. This should allow the formalism developed for the investigation of radiation-induced structural defects in nuclear reactor materials to be extended to the investigation of device performance degradation effects induced by particle radiation found in spacebased systems.

Fig. 1. Two-dimensional illustrations of different types of point defects in a crystal
The rest of the paper is organized as follows. In section II, we present our atomic-scale MD model to cover both the ultra-fast defect generation and intermediate defect stabilization stages, as well as the mesoscopic-scale rate theory for defect migration and interaction processes. In Section III, master equations for both planar dislocation-loop and spherical void growth are introduced for studying surface and bulk sink dynamics, respectively. In section IV, master equations are presented for exploring the steady-state spatial distribution of defects in layered structure materials. In addition, a density-functional theory is introduced for specifying electronic properties of point defects and four device physics models are employed for characterizing and understanding defect-assisted resonant tunneling, reduced carrier mobility, non-radiative recombination with defects and inelastic light scattering by charged defects. In section V, an indication of how this defect theory might be applied to experimental data describing device performance degradation is presented. Finally, some concluding remarks are presented in section VI.

Model and Theory
Atomic-Scale Modeling for Ultrafast Defect Generation (Displacement Cascade: t < 100 ps) A schematic of a displacement-cascade event by proton irradiation is shown in Fig. 2. For the neutronnucleus elastic collision, this process can be simply regarded as colliding hard spheres as an approximation due to their charge neutrality. The more complicated inelastic collision of neutrons with a nucleus, however, could involve generating an additional neutron [(n, 2n)process] or photon emission [(n, γ)-process], which are both important to the displacement of atoms. For the proton-nucleus elastic collision, on the other hand, the extra interaction (potential function) between the electron cloud and the proton should be considered.
Commonly, the end product of a particle collision results in the PKA having an excess kinetic energy and the subsequent atom-atom interaction represents the most fundamental physical mechanism of the radiation displacement damage (Gao and Weber, 2000a;Gao and Weber, 2003b). Since radiation damage events are random in nature, a large number of damage events are required to obtain good statistics by choosing different directions and locations for the PKA. On the other hand, the dynamics in the damage procedure can be accurately described by employing a realistic interatomic potential (MD model). For incident charged particles, the detailed form of the interatomic potential depends on the closest separation between two collision partners, which itself is determined by the kinetic energy of the incident particles (e.g., heavy-slow ions and relativistic electrons).
The point defect generation as a result of displacement cascades is closely related to the PKA energy, which can be described statistically by an average transfer energy to the PKA. Such an average transfer energy can be calculated by using the energyloss theory and measured by the so-called proton (electron) energy-loss spectroscopy (Gumbs, 1989) as a function of various incident charged particle energies. The defects can also be identified experimentally by using positron annihilation (Tuomisto and Makkonen, 2013). With help from the computed energy loss of incident particles per unit length (called the loss function), the range of the particle before its full stop inside a crystal can be found. On the other hand, the MD method has been widely employed to simulate defect generation in a number of semiconductors, including Si, (Rubia and Gilmer, 1995) SiC, (Gao and Weber, 2000) GaAs, (Nordlund et al., 2001) and GaN (Nord et al., 2003).

Fig. 2. Two-dimensional schematic of a displacement cascade induced by incident protons on a crystal
These simulations provide important insights into the mechanisms for defect generation in semiconductors and predict the number and type of defects, spatial distribution of defects and initial correlation among defect species produced by the incident radiation for subsequent device level models.
Basically, in MD simulations the time evolution of a set of interacting particles is tracked via the solution of Newton's equations of motion as shown below: where, the indices j = 1, 2,..., N label individual N particles in the system, r j (t) [x j (t), y j (t), z j (t)] is the position vector of the jth particle and ( ) is the force acting upon the jth particle at time t with interacting potential V jk between the jth and kth particle and m j is the mass of the corresponding particle. In general, F j (t) will depend on both particle positions and velocities at time t. Here, the energy loss of protons to electrons is ignored due to thin samples considered. To integrate the above second-order differential equations, the instantaneous forces acting on the particles and their initial positions and velocities need to be specified. Due to the many-body nature of the problem, the equations of motion have to be discretized and solved numerically. The MD trajectories are defined by both position vector r j (t) and velocity vector and they describe the time evolution of the system in position-velocity phase space. Accordingly, the positions and velocities are propagated with a small time interval ∆t using numerical integrators. The numerical integration of Newton's equations of motion is employed to find an expression that defines positions r j (t + ∆t) at time t + ∆t in terms of the already known positions r j (t) at time t. Because of its simplicity and stability, the Verlet algorithm is commonly used in MD simulations. (Frenkel and Smit, 2002) However, other popular algorithms, such as Leapfrog, Velocity Verlet, Beeman's algorithms, (Frenkel and Smit, 2002;Allen and Tildesley, 1987) predictor-corrector, (Gear, 1971) and symplectic integrators, (Tuckerman and Martyna, 2000) are also widely adopted. For non-PKA particles, their two initial conditions can be set as r j (−∆t) = r j (0) = R j , where R j is the lattice vector for the jth site. If the PKA is given an initial velocity v 0 , in addition to In MD simulations, the atomic force field is crucial to determine physical systems in which collections of atoms are kept together by interatomic forces that can be calculated from empirical or semi-empirical interatomic potentials. Because of extensive applications of MD methods in materials science, a variety of techniques have been utilized over the years to develop reliable atomic-potential models. One of the early successful attempts to include manybody effects was the introduction of the embedding functional, (Norskov and Lang, 1980) which depends nonlinearly upon the coordination number (defined below) of each atom. This development led to the birth of the Embedded Atom Method (EAM), (Finnis and Sinclair, 1984) which provides a relatively accurate description for noble transition metals as well as their alloys. However, the Tersoff potential formalism (Tersoff, 1988) is based on the concept of bond order and has been widely applied to a large number of semiconductors. Novel many-body forms have been tried in an attempt to capture as much as possible the physics and chemistry of the bonding. A typical analytical form is constituted by a number of functions, depending on geometrical quantities, such as distances or orientations, or on intermediate variables, such as atom coordinations. For example, a Tersoff potential has the appearance of a pair potential as below: where, the terms with i = j are excluded in the above summations, r ij = |r i -r j | and the first and second terms represent repulsive and attractive interactions, respectively. However, the second term in Equation (2) is not a true pair potential since B ij is not a constant. In fact, it is the bond order for the bond joining the ith and jth atoms and it is a decreasing function of a "coordination" G ij assigned to the bond. Therefore, we have B ij = B(G ij ) and G ij is in turn defined by: where, f c (r) and g(θ) are suitable functions. The basic idea is that the i-j bond is weakened by the presence of other i-k and j-k bonds involving the intermediate atom at site k. The amount of bond weakening is determined by where the other bonds are. Angular terms appear necessary to construct a realistic model. When using a potential, the researcher should always be familiar with its chemical transferability properties, such as bond length, bond angle and bond energy and validate critically the results obtained in unusual conditions, for example, for very low coordinations, very high temperature, or very high pressure.
Atomic-Scale Modeling for Intermediate Defect Stabilization (Stable Defect and Cluster Formations: 100 ps<t<10 ns) The atom-displacement-generated point defects (vacancies and interstitial atoms) under particle irradiation will thermally diffuse in space and interact and react with dynamically distributed bulk sinks, planar dislocation loops, (Hirth and Lothe, 1982) spherical voids and clusters (due to collision cascade) at the same time. Generally, the kinetic energy of the incident particles (or equivalently, the recoil energy of the struck atom) determines the species and number of individual point defects during the initial phase (in addition to the rate of defect generation), while the flux of the incident particles decides the defect density and the nature of point-defect diffusion, (Mehrer, 2007) i.e., either in an independent way (for lowdensity non-interacting point defects) or in a direction-correlated way (for highdensity interacting point defects) (Posselt et al., 2008). The macroscopic property changes of the irradiated system are related to the particle energy-flux per unit time by the so-called damage function, or factor, which is extracted by experimental measurements. However, the damage factor is found to depend on the initial approximation in a sensitive way. Therefore, we are not able to treat physically the radiation displacement damage effects as a black box through a fitting procedure. Instead, we should understand the full dynamics of these defects on all time scales after they have been produced. The spatial distribution of the mobile Frenkel pairs (i.e., vacancy-interstitial pairs) that are created is crucial in determining the number that survive annihilation or immobilization by clustering due to damage cascade.
The statistically-averaged spatial distribution of point defects that are generated can be calculated based on the defect formation and recombination rates, as well as the follow-up processes for defect diffusion, interactions and reactions (Posselt et al., 2005). If the degree of atom displacements is limited due to high incident particle kinetic energies and low number intensities, we generally seek the radiation degradation effects on electronic and optoelectronic devices rather than looking at material level radiation damage effects when there is a significant level of atom displacements under intense low-energy particle irradiation (Strour et al., 2003). This radiation degradation depends not only on the particle radiation source and material, but also on the device structure and functionality. The analytical theory below can only provide a qualitative understanding of the collision-and thermally-activated diffusion processes, while the MD calculation based on a realistic interatomic potential is able to provide a quantitative result for comparison with experimental data.

Point-Defect Generation Rate
The spatially-temporally-dependent damage rate per unit volume for the displacement atoms in a crystal can be calculated from : where, n at is the crystal atom volume density, I ext (t|ε i ) represents the external dynamical energy-dependent particle intensity per unit energy, σ D (r|ε i ) stands for both the position-and energy-dependent displacement cross section and E min (E max ) corresponds to the minimum (maximum) kinetic energy in the energy distribution of incident particles.
Since the displacement cross section σ D (r|ε i ) in Equation (4) physically describes the probability for the displacement of struck lattice atoms by incident particles, we can directly write down: where, σ C (r|ε i , ε R ) is the differential energy transfer cross section by collision, which measures the probability that an incident particle with kinetic energy ε i will transfer a recoil energy ε R to a struck lattice atom, N D (r|ε R ) represents the average number of displaced atoms due to collision and ε 1 (ε 2 ) labels the minimum (maximum) recoil energy acquired by the struck lattice atom. Although N D (r|ε R ) can be directly determined by MD simulation, (Terentyev et al., 2006) the simplest approximation to estimate the displaced atoms is the Kinchin-Pease model (for a solid composed of randomly arranged atoms by ignoring focusing and channeling effects). (Kinchin and Pease, 1955;Olander, 1976) In practice, N D (r|ε R ) can be accurately determined by either classic (Devanathan et al., 1998) or ab initio MD simulations (Gao et al., 2011). In addition, the magnitude of the differential energy transfer cross section σ C (r|ε i , ε R ) introduced in Equation (5), which can be regarded as the crystal response to an external particle collision with lattice atoms, depends on the detailed collision mechanism and the form of the scattering potential as well. Here, as a simple example, we first give the expression of the differential energy transfer cross section for the elastic scattering. For the well-known Rutherford elastic scattering model (Rutherford, 1911) based on an unscreened Coulomb potential U R (ρ) = Z 1 Z 2 e/ ∈ 0 ρ for protons with ρ being the radius in the local frame centered on the lattice atom, we get: where, γ(r) = 4mM(r)/[M(r) + m] 2 , m [M(r)] is the mass of the incident particle (surface atoms or different lattice atoms), b 0 (r) = Z 1 Z 2 (r) e 2 /η(r)∈ 0 ε i with Z 1 and Z 2 (r) being the nuclear charge numbers for particles and different lattice atoms and η(r) = m/[M(r) + m].
If the kinetic energy of the incident particles is very high, the Rutherford scattering model becomes no longer applicable. In this case, we have to consider hard-sphere type collisions for neutrons, which leads to: where, a Born-Mayer potential U B-M (ρ) = A exp(-ρ/B) is employed (Abrahamson, 1969). On the other hand, nuclear scattering with heavyslow ions represented by a power-law interacting potential U I (ρ) = (e/∈ 0 a B ) (Z 1 /Z 2 ) 5/6 (a B /ρ) 2 leads to: is the screening length with a B being the Bohr radius and E a (r) = (e 2 /∈ 0 a B ) [Z 1 /Z 2 (r)] 7/6 /η(r).
In the special case, for the incidence of relativistic light electrons, we have: where, β 0 = υ/c with υ being the velocity of incident , m 0 is the freeelectron mass and α 0 (r) = Z 2 (r)/137. Moreover, we have the relation ≤ for the relativisticparticle velocity and kinetic energy.
For the case of isotropic inelastic scattering with an energy loss Q 0 , we have the differential energy transfer cross section: where, A(r) = M(r)/m and σ is (r|ε i , Q 0 ) represents the isotropic differential energy transfer cross section for the resolved resonance in the center-of-mass frame.
In addition, the incident-particle kinetic energy is no longer a constant if the particles are charged, e.g., protons and ions. In this case, we have with ε 0 being the incident-particle energy at the left boundary z = 0, where a layered structure in the zdirection is assumed. As a result, we find that electronic stopping will dominate at short distances, while elastic collisions will dominate near the end of the range. As an example, by using σ C (r|ε i , ε R ) from Equation (6) and N D (r|ε R ) from the Kinchin-Pease model, the displacement cross section from Equation (5) becomes: where,  . Furthermore, by using the result in Equation (11), the displacement damage rate per unit volume from Equation (4) is: where, are the average values with respect to the incident particle intensity per unit energy L ext (t|ε i ) in the energy range of E th (r)/γ≤ε i ≤F ext (r,t) is the integrated external particle intensity for the same energy range and the term in the bracket is the number of Frenkel pairs produced per incident particle.

Point-Defect Diffusion Coefficient
Even in the absence of particle irradiation, there still exist some thermally activated vacancies at room temperature in a crystal. In this case, the Helmholtz freeenergy function in thermodynamics can be applied by assuming the volume of the crystal is a constant. In the presence of crystal defects, both the entropy S and the enthalpy H p of a perfect crystal will be changed. A straightforward calculation gives the thermal-equilibrium numbers of vacancies eq v C and interstitials eq i C as follows: where T is the temperature, E v is the vacancy formation energy, which is smaller than the interstitial formation energy E i and S v (S i ) is the change in entropy due to vibrational vacancy (interstitial) disorder. Diffusion of defects is driven by forces other than the concentration gradient of defects, such as stress or strain, electric fields, temperature, etc. The second Fick's law (Murty and Charit, 2013) directly gives rise to the following diffusion equation on the macroscopic scale: where, D v (r, t) = D(r, t|C v ) and D i (r, t) = D(r, t|C i ) are called the diffusion coefficients for vacancy and interstitial atoms, respectively and c v,i (r, t) = C v,i (r, t)/V is the defect concentrations with V being the volume of the system considered.
In addition, by assuming a microscopic random walk for the diffusion process, we get the Einstein formula (Guinan et al., 1977): where, the temperature-independent part, D 0 (r), is proportional to the Debye frequency (∼ 10THz) and is independent of defect concentration, E ac (r) is the activation energy for thermal diffusion λ d (r) is the diffusion length and Γ(r) is the defect jump rate. For tracer-atom (single radiative atom) diffusion, the random-walk model can not be used. Instead, the diffusion process becomes correlated, described by the Haven coefficient f(r) and we obtain where f(r)<1 depends on the crystal structure and the diffusion mechanism.
The lattice-atom correlated diffusion coefficients , ( , ) v i a D r t for vacancy and interstitial cases are given by: which depend on the defect concentrations in this case, implying a nonlinear diffusion equation. It should be noted that the correlation factors, f v,i (r) and diffusion coefficients, D v,i (r), can be determined by atomic-level simulations (Guinan et al., 1977;Posselt et al., 2008). Equation (16) can be directly applied to defect diffusion under irradiation, as long as the defect concentrations are known.

Mesoscopic-Scale Rate Theory for Slow Defect Migration and Interaction (Defect Reaction and Migration: t > 10 ns)
The formation, growth and dissolution of defect clusters such as voids, dislocation loops, etc., depend on the diffusion of point defects and their reaction with these defect clusters (Bullough et al., 1979;Colubov et al., 2000). At the same time, they also depend on the concentration of point defects in the crystal. Since particle irradiation greatly raises the defect concentration above its thermal-equilibrium value, the diffusion coefficient can be enhanced. It can also be enhanced by the creation of new defect species.

Point-Defect Diffusion Equation
By introducing the local coupling rates ( , ), ( , ) and Γ vs (r, t) for vacancy-interstitial recombination, interstitial-sink and vacancy-sink reaction rates, we can write down the following two nonlinear rate-based diffusion equations for binary crystals (Wiedersich et al., 1979): where, we have neglected correlated diffusions and defect-defect interactions. Moreover, J v (r, t|c v ) and J i (r, t|c i ) in Equations (17) and (18) are the particle currents for vacancies and interstitials and the master equation for determining the local sink concentration c s (r,i) will be given later. For simplicity, we consider a homogeneous system with volume V in the absence of vacancy and interstitial currents and assume all the rates are independent of time. We further neglect the small thermal-equilibrium vacancy concentration and write the defect generation rate as G 0 = G 0 V. For this model system, we find that the evolution of c v (t) and c i (t) depends on the temperature and c s and can be characterized in several regimes separated by different time scales τ, including initial buildup without reaction, dominant vacancy-interstitial mutual recombination and final vacancy and interstitial annihilation by sinks. Physically, it is easy to understand that the initial buildup and recombination regimes correspond to the 'ultra-fast' atomic-scale modeling, while the interstitial and vacancy annihilation regimes are associated with the 'slow' mesoscopic-scale modeling which may be simulated using the so-called rate theory (Singh and Zinkle, 1993) or phase field model (Li et al., 2010).

Recombination and Sink Annihilation Rates
In general, the reaction rate between species A and B can be expressed as Γ AB c A c B where c A and c B are the concentrations (particles/cm 3 ) and Γ AB (cm 3 /s) is the rate constant.
As an example, the recombination rate constant ( , ) r t ℜ in Equation (17) and (18) for vacancies and interstitials takes the form of (Waite, 1957): where, z iv (r) (an integer) is the bias factor, depending on the crystal structure and species, Ω(r) is the atomic volume, a 0 (r) is the lattice constant and D i (r, t) is the mobile interstitial diffusion coefficient. In a similar way, the interstitial-sink Γ is (r, t) or the vacancy-sink Γ vs (r, t) annihilation rate constants in Equations (17) and (18) as a a k r t c r t D r t = , where α corresponds to mobile defect species and κ αs (r, t) (cm −2 ) represents the sink strength given by (Doan and Martin, 2003): The sink strength measures the affinity of a sink for defects, which is independent of defect properties and corresponds to the mean distance for a traveling defect in the crystal before it is trapped by sinks.

Point-Defect Interaction Rates
In the absence of the macroscopic-scale gradient of defect concentration, the reaction between defects and sinks is reaction-rate-controlled. According to Equation (19), the defect-void interaction (Brailsford et al., 1976) can be described by the rate constants Γ {i,v}V (r, t) by: where, represents the radius of a void sphere involving n vacancies. The void strength is given by 2 2 0 where c V (r, t|n) is the concentration of voids containing n vacancies in the crystal.
Similarly, for the defect-dislocation interaction, we have the rate constants Γ {i,v}d (r,t) (Heald and Speight, 1975) (in units of cm 2 /s) given by: where, we replace Ω in Equation (19) by an atomic area 2 0 ( ), ( ) ( ) id vd a z r z r ≠ ∼ and ρ dr,t (r, t) is the dislocation areal density.
Reactions driven by defect concentration gradients are diffusion limited instead of reaction-rate limited as discussed above. In this case, we have to solve the diffusion term ∇· [D i,v (r, t) ∇ ci,v (r, t)] with the generation term G 0 (r, t) for spherical (voids) or cylindrical (dislocation lines) coordinates. In the case of dislocation loops, the sink strength is more complicated and increases with increasing loop size (Bullough et al., 1979).
For the defect-void interaction, we get . For the defect-dislocation line interaction, on the other hand, we have where R 0 is the absorption radius of a sink, R id (r, t) and R vd (r, t) are the sink capture radii for interstitials and vacancies, respectively with id vd R R ≫ .
In the presence of Grain Boundaries (GBs), we obtain the grain-boundary sink strength: where, R gb (r) is the radius of a spherical grain, c gb (r) is the grain concentration and κ i,v (r, t) is the sink strength for the grain interior for interstitlals or vacancies due to dislocations and voids. Moreover, its rate constant is More recently, atomic-level simulations have been employed to study the sink strength of grain boundaries (Tschopp et al., 2012;Jiang et al., 2014). These simulations revealed that the sink strength can exceed the theoretical strength of GBs due to the effect of GB stress fields.

Radiation-Induced Segregation
For a binary A-B alloy (or donor and acceptor randomly-doped semiconductors), in the absence of sinks, the diffusion equations for vacancies, interstitials and atoms A and B are (Wiedersich et al., 1979): where, are the diffusivity coefficients and the dimensionless X(r,t) is the thermodynamic factor connecting the chemicalpotential gradient to the concentration gradient (or built-in electric field at junctions). In addition, we have c B (r, t) = Ω 1 (r) -c A (r, t) when small defect concentrations are neglected.
By requiring J A = J B = 0 and J i = J v for steady state and neglecting G 0 (r, t) and ( ) , r t ℜ in Equation (24)-(26), we get: where, N A,B = c A,B Ω and the direction of ∇c A can be either parallel or anti-parallel to ∇c V . Additionally, the undersized (oversized) solutes bounded to interstitials will be concentrated (depleted) around sinks to create a concentration gradient after their redistribution.
On the other hand, the oversized or undersized solutes with respect to the lattice atoms can act as traps for vacancies and interstitials, including release of defects from traps, recombination with trapped defects, trapping of free point defects and loss to internal sinks. This is further supplemented by three rate equations for trap and trapped defect concentrations.

Sink Dynamics
The growth of dislocation loops and spherical voids is determined by solving the point-defect balance equations without diffusion terms. Since the defect concentration is still changing with time due to the timedependent radiation source (or defect production rate), only quasi-steady state can be defined for short periods of time. Physically, the quasi-steady state is related to the fact that the change in sink strength due to microstructure evolution is slow compared to the response time of the defect population.

Planar Biased Dislocation-Loop Growth
By defining the dislocation line direction s and the Burgers vector b for edge (b ⊥ s) or skew (b||s) dislocations, the Peach-Koehler equation (Weertman, 1965) gives us the force f per length as: where, σ is the stress tensor. The force f along the b direction is the glide force, while f perpendicular to both s and b directions is called the climb force. The Peach-Koehler equation can be used for calculating the interaction between dislocations, where b and s are assigned to the target dislocation while σ is for the source dislocation. For an edge dislocation, we have five non-zero stress tensor elements σ xx , σ yy , σ zz and σ xy = σ yc , while for a skew dislocation, we only have four non-zero stress tensor elements σ yx = σ zy and σ xz = σ zx . Besides the dislocation lines, there also exists Frank loops. For example, the close-packed fcc lattice follows the stacking sequence ABCABCABC …, where A, B and C correspond to different planes of atoms. It can be modified to ABCAB/ABC …, where "/" denotes the intrinsic single fault or missing plane of atoms. It can also be modified to ABCAB/A/CABC …, where a plane of atoms or the extrinsic double fault is inserted.
Interstitial condensation can occur around the edges of the depleted zone. A cluster of point defects can be a line, a disc or a void. The formation of perfect or faulted loops of interstitials competes with the formation of voids of vacancies, which are also affected by the irradiation temperature.
The nucleation of loops is a clustering process that results in a critical size embryo for further growth. As an example, by denoting the number of clusters consisting of j vacancies as ρ v (j), the master equations for ρ v (j) are: where, β vn and β in are the capture rates of migrating vacancy (v) or interstitial (i) clusters of size n by a cluster of size j and ( ) n v a j is the emission rate for the new vacancy cluster of size n by a cluster of size j. In Equation (29), the first term is the direct production of a cluster of size j, while the second term is the loss of clusters from size j due to absorption of a cluster of size n. The third term is the loss of a cluster of size j due to emission of a cluster of size n. The fourth and fifth terms in Equation (29) are the addition of clusters to the cluster of size j due to absorption of vacancy clusters by a smaller cluster and absorption of interstitial clusters by a larger cluster and the last term is the addition of clusters to the cluster of size j due to loss of vacancy clusters by a larger cluster.
Since the dominant contribution for cluster reactions is with point defects (i.e., cluster of size j = 1), for both vacancies and interstitials, Equation (29) with j≥2 can be simplified to: If the cluster size index j can be treated as a continuous variable ξ, Equation (30) reduces to a Fokker-Planck equation (Semenov and Woo, 2002).
According to Equation (30), for the dislocation loop growth, we find the evolution of the number density ρ il (j, t) for the interstitial loop of size j satisfies: Where: r(j) and z c are the radius and bias factor of an interstitial loop of size j and E b (j) is the binding energy for a cluster of j defects.
The saturation of the dislocation density ρ d (t) in quasi-steady state was found experimentally to be due to a recovery process at high temperatures,  given by: where, A 2 and A 1 are constants. This gives the steadystate solution and the time-dependent solution is found to be: where, 0 d ρ is the initial value and In addition, a phase field model has been developed to describe the growth kinetics of interstitial loops in irradiated materials during aging. The diffusion of vacancies and interstitials and the elastic interaction between interstitial loops and point defects are accounted for in this model. The effects of interstitial concentration, chemical potential and elastic interaction on the growth kinetics and stability of interstitial loops are investigated in two and three dimensions. The elastic interaction enhances the growth kinetics of interstitial loops. The elastic interaction also affects the stability of a small interstitial loop adjacent to a larger loop. The linear growth rates for interstitial loops predicted is in agreement with the previous theoretical predictions and experimental observations (Hu et al., 2012).

Spherical Neutral Void Growth
Not all the defects generated by radiation-induced atom displacement are point defects. Some of the defects form clusters and the vacancy clusters may further grow to form voids. For a small number of vacancies, the spherical void is favorable, while for a large number of vacancies, the planar loop is a more stable configuration.
The dynamics for void growth is very similar to that for dislocation-loop growth. The net absorption rate of vacancies by a spherical void is the difference of absorption rates of vacancies and interstitials, i.e., Therefore, the equation for the growth of a spherical void of radius R V (t) (or volume) in quasi-steady state is: where, c v (R V ) is the vacancy concentration at the void surface.
From the balance equation, we get the concentrations of point vacancies and interstitials as follows: where, Inserting Equation (37) into Equation (36), we obtain: where, dR 0 (t)/dt is proportional to (z i − z v ) G 0 and independent of temperature, while the second negative term represents the thermal emission of defects from sinks and strongly depends on temperature (proportional to eq v v D C ). Neglecting the thermal emission in Equation (38), we get from Equation (36): where, the sign of the bias of dislocation ( for vacancies and interstitials determines the occurrence of either growth (dR V /dt > 0) or shrinkage (dR V /dt < 0). It should be noted that void growth in materials under radiation has been recently studied by the phase field model, (Li et al., 2010) which provides important insights into the growth kinetics of voids. The model takes into account the generation of vacancies and interstitials associated with the irradiation damage, the recombination between vacancies and interstitials, defect diffusion and defect sinks. The results demonstrate that the temperature gradient causes void migration and defect fluxes, i.e., the Soret effect, which affects void stability and growth kinetics. It is found that the effect of defect concentration, generation rate and recombination rate on void mobility for migration is minor although they strongly influence the void growth kinetics.

Radiation Degradation of Electronic Devices
Let us consider a commonly used layered-structure material, (Huang and Zhou, 1988) as shown in Fig. 3. Each material layer is characterized by the radiation parameters G j , j ℜ , D j and Γ j with j = 1, 2, 3, 4 for generation and recombination rates, diffusion coefficient and bulk-sink annihilation, which will be employed to model the dynamics from an ultra-fast atomic-scale up to 100 ns. The calculated non-steady state defect distribution in each layer will be used for initial conditions in a slow mesoscopic-scale diffusion and annihilation model in order to calculate the steady-state spatial distribution of defects in the whole layered structure. In modeling the mesoscopic-scale, the interface-sink strengths 2 i k with i = 1, 2, 3 will also be considered. Once the steadystate distribution of point defects, denoted as ρ d (z), is obtained for the whole layered structure, they will be fed into the follow-up calculations for radiation degradation in electronic devices, as described below.
The band structure of a crystal largely determines the properties of electrons, (Gumbs and Huang, 2013) such as effective mass, bandgap energy, density of states, plasma frequency and absorption coefficient. These electron properties are a result of the unique crystal potential from all lattice atoms, instead of properties of an individual lattice atom. On the other hand, the radiation-induced displacements of lattice atoms are determined not only by the intrinsic properties, such as mass of the atoms, but also by the extrinsic conditions, such as kinetic energy of incident particles and lattice temperature.

Steady-State Defect Distributions
For the reaction rate control system shown in Fig. 3, by generalizing Equations (17) and (18), we write down the diffusion equations for point vacancies and interstitials as:  Fig. 3. Layered-structure material with radiation parameters G j , R j , D j and Γ j (j = 1, 2, 3, 4) for generation, recombination, diffusion coefficient and bulk-sink annihilation, respectively. In addition, 2 i k for i = 1, 2, 3 represents the interface-sink strengths. Particles are incident from the front surface at z = 0 and exit from the back surface at z = z 4 where, the integer j is the layer index, z j-1 and z j represent the left and right interface positions of the jth layer, respectively. In Equations (40) and (41), we used the facts that for a reaction rate control system.

The diffusion coefficients
where, , 1 j i v f < is the diffusion correlation factor, η j is the structural factor relating to the jump distance and number of nearest neighbors, ν j is the jump frequency on the order of the Debye frequency,  (40) and (41) is determined from the following void growth equation: c t R is the vacancy concentration at the void surface and is given by: γ j is the surface tension of the void, σ j is the hydrostatic stress applied to the void and 0 ( ) v c j is the thermalequilibrium vacancy concentration, given by: is the change of entropy for the formation of a point vacancy and j f E is the point vacancy formation energy. By using a continuous variable, the void concentration ( , ) j V c z t n ( 2) with n ≥ introduced in Equations (40) and (41) can be obtained by solving the Fokker-Planck equation in the size space below (with ξ = n): where, is the cluster production rate per volume (with ζ ≥ 2). Finally, the dislocation-loop density ( ) j d t n ρ introduced in Equations (40) and (41) can be found from (with n ≥ 4): a z t n t n z t n t n z t n z t n a z t n t n β ρ ρ β ρ β β ρ where, G 0 (z j , t|n) is the production rate per unit volume for the interstitial dislocation loop of length n, the absorption (β i,v ) and emission (α i,v ) rates in Equation (47) are defined by: E n is the binding energy for a cluster of n interstitial atoms.
The initial condition for the diffusion equations will be given by the corresponding calculated results from the atomic-scale model for individual layers. The point defect diffusion occurs mainly around interfaces between two adjacent layers or across the interfaces. The boundary conditions with continuous concentrations of point defects, as well as the jump in their derivatives determined by the dislocation sinks, will be applied at each interface. In addition, the constraints for the zero concentration of point defects as well as the zero derivative of the concentration with respect to z at the two surfaces of the system will be enforced in our numerical computations.

Point-Defect Electronic States
To study the defect degradation effect on devices, we need to know not only the concentration and spatial distribution ρ d (r) of the irradiation-induced defects but also their electronic properties, such as energy level E j , wave function ψ j (r) and local density of states D d (r, E). Although the semi-classical MD calculation and the reaction-rate theory allow us to obtain the concentration and spatial distribution of defects, we still require Density-Functional Theory (Drabold and Estreicher, 2007;Freysoldt et al., 2014) (DFT) for calculating defect configurations, energy levels, density of states and charge trapping by point defects in crystals.
The main idea of DFT is to reformulate the energy of an atomic system as a functional of the ground state electron density function ρ 0 (r) instead of individual electron wave functions. The proof of existence of such a functional relies on a one-to-one correspondence between the external potential V ext ({R ℓ }, {r m }) and ρ 0 (r), where {R ℓ } and {r m } label all the lattice atoms and electrons, respectively. The mapping of V ext ({R ℓ }, {r m }) onto ρ 0 (r) is obvious. Any Hamiltonian H ⌢ with a given external potential V ext ({R ℓ }, {r m }) has a ground state solution with an N-electron wave function ϕ 0 ({r m }), which can be uniquely associated with the electron density function ρ 0 (r) using: Due to the resulting one-to-one correspondence between V ext ({R ℓ }, {r m }) and ρ 0 (r), the energy E i of the atomic system can be expressed as a functional of the electron density ρ 0 (r). The many-electron wave function of ϕ(r 1 , r 2 ,…, r N ) depends on the 'combination' of all spatial electron coordinates. Unfortunately, such an approach would by far exceed any computational capabilities. However, this problem can be overcome by using the Kohn-Sham (KS) ansatz, (Kohn and Sham, 1965) in which the fully-interacting system is replaced by a non-interacting one. This approach corresponds to a mean-field approach, where the many-electron wave function is decomposed into a product of N singleelectron orbitals φ i (r) (i.e., Slater determinant). This simplification leads to a neglect of an energy contribution termed 'correlations'. As a correction, the functional E xc [ρ(r)] must be introduced as an additional term in the Hamiltonian. Applying the variation principle to the modified Hamiltonian yields a single-particle-like Schrodinger equation, also referred to as the Kohn-Sham equation in DFT. This equation includes an effective potential V eff (r), which is produced by the Coulomb forces of all other electrons and nuclei and incorporates the exchange and correlation interactions, i.e.,: where, V ee (r) describes the electron-electron interaction (the classical Coulomb interaction) that is defined by: and V xc [ρ(r)] is the functional derivative of the exchange correlation energy with respect to the electron density function: The total energy of the atomic system can be arranged as: where, T k [ρ] represents the kinetic energy of noninteracting electrons. The exchange-correlation functional in Equation (54) can be written as: where, V el [ρ] is non-local electron-electron interaction beyond the classical one in Equation (53). E xc [ρ] is simply the sum of the error made in using a noninteracting kinetic energy and the error made in treating the electron-electron interaction classically. The Kohn-Sham equations in Equation (51) have the same structure as the Hartree-Fock equations with the nonlocal exchange potential replaced by the local exchange-correlation potential V xc [ρ(r)]. The computational cost of solving the Kohn-Sham equations scales formally as N 3 (due to the need to maintain the orthogonality of N orbitals), but in current practice it drops to N through the exploitation of the locality of the orbitals. Actually, the utility of the theory rests on the approximation made for E xc [ρ].
Therefore, the correct description of the exchangecorrelation functional takes a crucial role in DFT. The localdensity approximation has already achieved satisfactory results for systems with a slowly varying electron density function, such as metals (Kohanoff, 2006). However, it has a tendency, termed over-binding, to overestimate binding energies and thus, for instance, predicts too strong hydrogen bonds with too short bond lengths. The generalized gradient approximation is a systematic expansion, which gives good results in most cases and corrects the issue of over-binding (Kohanoff, 2006). Recently, hybrid functionals  have emerged, which achieve an improved accuracy, especially for semiconductors with a bandgap .
Defect levels for charge capture or emission are calculated by means of the formation energies ' [ ] q q f E X , (de Walle and Neugebauer, 2004) with are defined for a certain charge state q and a certain atomic configuration ' q X of the defect as: Here, E tot [bulk] stands for the total energy of a supercell containing pure bulk material while represents the super-cell containing a defect. The third term corrects for the different numbers of atoms in both super-cells. The integer n j stands for the number of added (n j > 0) or removed (n j <0) atoms which are required to create the defect from a perfect bulk structure and ξ j denotes the corresponding energy in an atomic reservoir, which must be specified for each individual case. The fourth term accounts for the charge state q of the defect, in which µ is defined as the electron chemical potential referenced with respect to the valence band edge ε ν in a bulk-like region and ∆V corrects the shift in the reference level between charged and uncharged super-cells and is obtained from the difference in the electrostatic potential far distant from the defect for these two super-cells. Due to the periodic boundary conditions, charge neutrality must be maintained within a super-cell. Therefore, a homogeneous compensating background charge must be introduced in calculations of charged defects. This artificial Coulomb interaction is corrected by the last term E corr .

Defect-Assisted Resonant Tunneling
At low temperatures, the defect-assisted tunneling through thermal emission can be neglected . Therefore, the whole elastic tunneling process can be divided into two subsequent ones, i.e., tunnel capture and tunnel emission, as shown in Fig. 4. Although the in-plane momentum of electrons is not conserved during the tunneling process, the kinetic energy of electrons is conserved. For a neutral point defect, let us assume that it sits at an arbitrary position z = z 0 inside a barrier layer (0≤z 0 ≤L B ) between the Left (L) and Right (R) electrodes with energy levels 0 < E d (z 0 ) <∆E c , where ∆E c is the conduction band offset for the middle barrier layer. A bias field ε b is applied across the layer, leading to a voltage drop Fermi distribution f L is first captured from left (process-1 in red) by the point defect through tunneling and then emitted to a continuum state above the energy barrier on the right (process-2 in blue) through tunneling in the presence of a voltage drop V b across the barrier layer By assuming a large voltage drop, we need to consider only the forward current from left to right but not the backward current from right to left. In this picture, the left-going (capture) tunneling current density J L (V b , T ) can be formally written as : where, ρ d (z 0 ) represents the distribution of the pointdefect concentration, U d (r) is the Coulomb potential associated with the point defect, ψ d (r) is its wave function and Γ d is the broadening in the density of states for point defects.
In Equation (58), the occupation factor F L (E k ) is defined as: where, E k = h 2 k 2 /2m * is the electron kinetic energy with effective mass m * in the left electrode, Z e and Z f represent the structural degeneracy factors of the point defect when empty or filled, g[E d (z 0 )] is the defect occupancy function and is the Fermi distribution function in the left electrode with chemical potential µ 0 . In addition, by employing the WKB approximation (Gill, 1982) for the electron wave function Ψ k (r), the interaction matrix 〈Ψ k |U d |ψ d 〉 is calculated as (Stievenard et al., 1992): where, A k is an unknown coefficient to be determined by the continuity of the wave function at the boundaries, S is the cross-sectional area: In a similar way, we can also calculate the rightgoing (escape) tunneling current density J R (V b , T ). In steady state, we have This allows us to eliminate the unknown defect occupancy function g[E d (z 0 )] and eventually obtain (Stievenard et al., 1992): where, . In addition, the tunnel-capture rate (probability) P c (z 0 ) of an electron by a point defect in Equation (63) is defined as: and the tunnel-emission rate (probability) P em (z 0 ) of an electron captured by a strongly-localized point defect is given by: For photo-detectors, the defect-assisted resonant tunneling greatly increases the dark current in the absence of incident light, which generates excess noise and reduces the detectivity of the photo-detector (Huang et al., 2000).

Reduced Carrier Mobility
When point defects are charged with a charge number |Z * | ≥ 1, they can scatter conduction electrons through their Coulomb potential  where the defect is charged and has an effective charge number Z * . The incident electron with wave vector k and kinetic energy E k is scattered into a different direction with wave vector k′ and kinetic energy E k′ . The scattering angle between k and k′ is denoted by θ kk′ and the elastic scattering process requires E k = E k′ In this case, the interaction matrix is calculated as (Huang et al., 2005a): where, the two-dimensional Fourier transform of U c (r-r i ) is denoted as Uc(q || , z-z i ) and give by: and ∈ r is the dielectric constant of the quantum-well host material. Moreover, q s in Equation (67) is the inverse Thomas-Fermi screening length for quantumwell electrons, given by (Huang and Manasreh, 1996): where, T is the electron temperature and µ 0 is the chemical potential of electrons in the quantum well.
Since the positions of point defects are random, by introducing a continuous linear density distribution ρ 1d (z 0 ) = S ρd (z 0 ) for point defects, the interaction matrix from Equation (66) becomes: where, L W is the width of the quantum well. Once the scattering matrix elements in Equation (70) are computed, by using Fermi's golden rule, the momentumrelaxation time τ 0 can be obtained from: where, ' || || k k θ represents the angle between the two inplane scattering wave vectors ' k and k . By using the momentum-relaxation τ 0 in Equation (71) The reduced mobility of conduction carriers by radiation-induced point defects will directly affect the speed of high-mobility field-effect transistors in an integrated circuit (Ando et al., 1982).

Non-Radiative Recombination with Defects
After the electrons are photo-excited from valence band to conduction band in a semiconductor, some of these photo-electrons will be quickly captured by point defects through an inelastic scattering process,  as shown in Fig. 6. By including the multi-phonon emission at room temperature (Huang et al., 2005b;Ridley, 1978), in this case the capture rate is calculated as (Jim´enez-Molinos et al., 2001;Garetto et al., 2012) where, U {e, h}p (r) represent the potentials for the electron-phonon and hole-phonon coupling, ψ d (r) is the wave function of the point defect, Γ d is the level broadening of the defect state,   in Equation (72) can be evaluated by (Huang and Alsing, 2008): where, V is the system volume, ∈ ∞ and ∈ s are the highfrequency and static dielectric constants of the host semiconductor and the inverse Thomas-Fermi screening length Q e,h for bulk electrons and holes is given by: where, ρ 1d (z 0 ) = Sρ d (z 0 ) is the linear density distribution of point defects in a layer with the crosssectional area S, g[E d (z 0 )] in Equation (78) The change in the non-radiative time by point defects in the system will reduce the quantum efficiency of photo-excited electrons in both light-emitting diodes and photo-detectors (Huang and Lyo, 1999).

Inelastic Light Scattering by Charged Defects
Let us choose the z direction perpendicular to the layered material. Light is incident on the layers in the xyplane and scattered by charged point defects within the layers. We consider incident light with photon energy hω i and wave vector k i scattered inelastically by bound electrons within point defects at ( , ) j j j r r z = for j = 1,2,….
If the scattered-light photon energy and wave vector are denoted by hω f and k f , respectively, the excitation energy and momentum transfer to charged point defects are given by  (Zhang et al., 1991): , e i and e f are the unit polarization vectors for the incident and scattered light, N ph (w) = [exp(hw/k B T) -1] −1 is the photon distribution function, T is the temperature, Ω q represents the solid angle in threedimensional q-space and ρ 1d (z 0 ) is the linear density of charged point defects.
In addition, the form factor ' 0 ( , ) nn z semiconductor) Read-Out Integrated Circuit (ROIC) and the detector array. FPAs developed for imaging applications in a space environment are meticulously characterized to evaluate their sensitivity, uniformity, operability and radiation hardness. Infrared (IR) hybrid detector arrays operated in the space environment may be subjected to a variety of radiation sources while in orbit. This means IR detectors intended for applications such as space-based surveillance or space-situational awareness must not only have high performance (high quantum efficiency and low dark current density), but also their radiation tolerance, or ability to withstand the effects of the radiation they would expect to encounter in a given orbit, must also be characterized. Pursuit of these performance metrics are intended to enhance technology development in multispectral sensors, higher pixel sensitivities, larger array formats with higher pixel densities, higher operating temperatures, simplified manufacturing techniques and reduced costs. Specific to space applications, as described above, is radiation hardness or tolerance, which is required to mitigate detrimental effects on IR detectors exposed to energetic charged particles (mostly electrons and protons and some heavy ions) confined by the earth's magnetic field (the van Allen radiation belts), galactic cosmic ray particles (which can be a range of elements or electromagnetic radiation) and particles, protons and heavy ions from solar flares. Detector technologies that operate in the harsh radiation environment of space with better radiation tolerance afford greater flexibility in orbit selection, technical applications and system sustainability and thus, are of more value to the space-based sensing community.
While specific mission requirements for rad-hardness vary, the effects of proton interactions with hybrid detector arrays tend to dominate in space. Thus, a specific detector's degree of radiation tolerance is typically first characterized by measuring its performance degradation as a function of proton irradiation fluence (Hubbs et al., 2007). Subjecting IR detectors to protons will lead to two distinct damage mechanisms, ionization damage and displacement damage, which are expected to account for the performance degradation due to irradiation that will occur in optical sensors in satellites. Ionization damage, or so-called Total Ionizing Dose (TID), effects occur as the incoming particles incident on the FPA give up their kinetic energy to ionizing additional electron-hole pairs within both the detector material and the CMOS silicon of the ROIC. Some of these excess carriers can become trapped in surface states and defect levels of the dielectric materials used for passivation and as gate oxide layers. These surface states and defect levels are to be determined with the theoretical approach presented in this study. The excess trapped charge manifests in flatband voltage shifts of the ROIC's MOSFETs (metaloxide-semiconductor field-effect transistors) and excess leakage currents in the detector pixels, typically surface-currents. Displacement damage results when the incoming particle's energy is lost to elastic or inelastic Coulomb scattering with an atomic nucleus that is sufficient to knock the atom from its lattice site and generate a vacancy-interstitial pair. The defect complexes so formed can act as electrically active sites where electron-hole recombination might occur (Claeys and Simoen, 2013). The additional defect generation thus shortens the detector material's minority carrier recombination lifetime, resulting in increased dark current, decreased responsivity and overall degraded uniformity.
Radiation tolerance characterization typically includes determining the rate of performance degradation via a damage factor analysis (Cowan et al., 2012;Morath et al., 2015). The degradation rates of each measured parameter of the detector's performance [dark current-including the diffusion-limited, Shockley-Reed-Hall (SRH) generation-recombination and trap-assisted tunneling mechanisms; lateral optical collection length or effective diffusion length-and thus minority-carrier recombination lifetime and mobility; quantum efficiency; etc.] are determined by plotting them as a function of the proton irradiation fluence Φ P and characterizing the changes. When the change in parameter X appears roughly linear with Φ P , which may be true on average or for only a certain fluence range, then a damage factor, K X , can be defined such that X = X 0 ±K X Φ P , where X 0 is the un-irradiated value of the performance parameter and the ± is determined by the expected change, e.g., + for dark current (Fig. 7) and −for quantum efficiency (Fig. 8). Damage factors are assumed to be related only to changes due to the effects of displacement damage, not ionization damage and are dependent on the proton energy E (Claeys and Simoen, 2013).
The radiation-induced defects can manifest in lower quantum efficiency (η) and higher dark current density [J D which is related to the tunneling current density in Equation (63)], due to the consequent reduction in minority carrier recombination lifetime [τ R which is related to the nonradiative recombination time τ nr in Equation (77)] and increase in SRH generation-recombination, carrier diffusion and trap-assisted tunneling. Thus, a substantial reduction in overall detector sensitivity, or "detectivity" (D*), is expected. Radiation tolerance can be characterized from performance measurements taken versus fluence by calculating the above-described damage factor (K X ) or the rate of degradation for each performance metric X (e.g., X = η, 1/τ R , J D , D*, etc). For comparison purposes, it is worth noting that damage factors are specific to the particle type and energy of the incoming radiation. With a known energydependence, K X (E), predictions of the expected onorbit degradation ∆X ideally become possible, according to: where, dΦ P (E)/dE is the expected orbit's differential proton fluence spectrum (Hubbs et al., 2007). Comparing damage factors from different detectors may allow for the more rad-tolerant detectors to be determined (as long as several other conditions that may impact damage formation are kept constant). A connection between the damage factor and the amount of displacement damage can be established by first assuming a linear increase in the defect concentration N T [related to the point vacancy concentration c v introduced in Equation (40)] with Φ P (which the data in (Hubbs et al., 2007;Cowan et al., 2012;Morath et al., 2015) support for the fluence ranges of interest). When dominated by SRH recombination, the minority carrier recombination lifetime τ R , assuming a single active recombination level, is given by: where, σ is the capture cross-section and ν th is the carrier thermal velocity. With these assumptions in mind, the minority carrier recombination lifetime must depend on Φ P according to: where, 1/τ R0 = N T (0) σ ν th . Based on the assumptions above and from Equations (90) and (91), a material's defect introduction rate is thus given by: Which means that if σν th is fairly constant, as expected for the fluence ranges of interest, K 1/τR is a relative measure of dN T /dΦ P , which is a true measure of a device's defect tolerance. Thus, Equation (92) suggests that if the various performance damage factors (for diffusion-limited or SRH or trap-assisted tunneling dark current densities, lateral optical collection length, quantum efficiency, etc.), can be directly related to K 1/τR , then an understanding of how they might relate to dN T /dΦ P and, furthermore, to each other, can be determined (Morath et al., 2015). The theory presented in this study will be used to determine the formation and dynamics of the proton-induced defect concentration and defect introduction rate that produce the damage factors that describe the degradation of sensor performance parameters. The experimentally observed increase of photo-detector dark current and decrease of quantum efficiency can be explained by defect-assisted tunneling and defect-induced capture models in Section IV, respectively.

Conclusion
In conclusion, we have proposed a multi-timescale microscopic model for fully characterizing the performance degradation of electronic and optoelectronic devices. In order to reach this goal, we will employ realistic interatomic potentials in a molecular-dynamics simulation for both the ultra-fast displacement cascade stage and the intermediate defect stabilization and cluster formation stage. This simulation will then be combined with a rate-diffusion theory for the slow defect reaction and migration stage. Additionally, with assistance from a density-functional theory for identifying defect species and their electronic properties, the calculated steady-state spatial distributions of defects and clusters will be used to study and understand the physical mechanisms that degrade electronic and optoelectronic devices, including defect-assisted resonant tunneling, reduced carrier mobility, non-radiative recombination with defects and inelastic light scattering by charged defects.
In this study, we have discussed several techniques for defect characterization. However, there are many other approaches for characterizing defect effects. These include electrical characterization techniques, such as deep-level transient spectroscopy and capacitancevoltage profiling and optical characterization techniques, such as cathodeluminescence and reflectance modulation. Physical and chemical characterization techniques can also be applied, including electron energy loss spectroscopy, (Gumbs and Horing, 1991) secondary ion mass spectrometry (Zhang et al., 2014) and chemical milling.
The molecular dynamics model presented in this study can be combined with a space-weather forecast model (Moldwin, 2008;Cooke and Katz, 1988) which predicts spatial-temporal fluxes and particle velocity distributions. With this combination of theories, the predicted irradiation conditions for particular satellite orbits allow electronic and optoelectronic devices to be specifically designed for operation in space with radiation-hardening techniques (such as self-healing and mitigation). This approach will effectively extend the lifetime of satellite onboard electronic and optoelectronic devices in nonbenign orbits and greatly reduce the cost. In addition, by improving the physical model for scintillation detectors, the accuracy of space-weather measurements will be enhanced. the other authors have read and approved the manuscript and no ethical issues involved.