Model for transport and reaction of defects and carriers within displacement cascades in gallium arsenide

Irradiation with energetic ions or neutrons causes atomic displacement cascades of lattice defects. In semiconductors, these defects capture electrons and holes. This increases the rate of carrier recombination and degrades the performance of minority-carrier devices such as heterojunction bipolar transistors (HBTs). This performance degradation may evolve with time as the defects diffuse and react. The objective of this work is to develop an atomistic description of carrier recombination at such displacement damage that includes the physics of the defects and their interactions with charge carriers, dopants, and each other. The model described here thus provides a microscopic description of the transport and reaction of charge carriers and defects within displacement cascades. The principal reactions of mobile defects and charge carriers are treated, in the continuum approximation, for a statistically representative defect cluster. The ultimate goal is to use this atomistic carrier recombination model as a sink term in macroscopic models of carrier flow within devices such as HBTs to simulate transient device response to damage.

Here, we describe an implementation of this recombination model for gallium arsenide (GaAs). A similar model was previously developed for silicon.¹ Such an atomistic model is more challenging for GaAs than for silicon, because the defect physics is more complex and less well understood. However, recent theoretical work using density-functional theory (DFT) has greatly improved knowledge of defects in GaAs. This new information on defect properties provides a basis for extending the atomistic model to GaAs. Evolution of defects and carrier recombination is considered for the example of p-type carbon-doped GaAs.

This paper is organized as follows. Section II describes the formation of defects in displacement cascades from energetic recoils and describes methods used to characterize defect clustering. Section II also describes the defect species included in the model and gives values for their properties, including formation and migration energies which were obtained from DFT simulations. Section II also gives the equations for transport and reaction of carriers and mobile defects on which the model is based. Section III describes computational methods used to solve the equations. Section IV presents results from simulations for conditions relevant to carrier recombination in III-V HBTs.

Energetic recoils from neutron collision events produce atomic displacement cascades with defects spatially localized in clusters. This clustering of the defects has profound effects on carrier recombination and defect annealing and is therefore a central feature of the model. The high local defect concentrations increase the rate of defect reactions and annealing. Furthermore, defects within a cluster may become charged producing a high electric field within the cluster which increases the flow of carriers into the cluster. Although neutron or ion irradiation produces collision cascades that vary greatly in size and structure, the clustering is approximated in the model as a single, representative, radially symmetric cluster.

The radial distribution of defects within the defect cluster is an initial condition in the simulation, and is determined from a pair correlation function analysis¹ of defect maps calculated using the binary collision code MARLOWE.² Such a correlation-function approach is fully objective, and avoids subjective choices relating to the definitions of a cluster and its center. The pair-correlation function is defined as

where ab represents pairings between defects of types a and b, c_i is the concentration of species i expressed as number per unit volume, and $r→$ is the spatial coordinate. In accord with the approximation of spherical symmetry, we average over orientation

where dA is the area differential on the surface of the sphere. In order to facilitate extraction of the correlation function from discrete-defect maps, we define

which is the average number of pairs with separation less than R per defect of type a. Inversion of Eq. (3) gives the finite-difference relation

Evaluation of the average radial concentration profile from the defect maps proceeds as follows. A separation distance R is calculated for each pair of defects in each cascade for a large number of recoil events. The function N_ab(R) is calculated from this list of pair separations and the pair correlation function $σ¯ab(R)$ is then evaluated using Eq. (4). For simulations discussed below, the MARLOWE calculations used a distribution of recoil energies for neutrons with an energy spectrum for a fast-burst reactor,³ as illustrated in Figure 1. The recoil energy distribution, also shown in Figure 1, was calculated using the Empire nuclear reaction code.⁴ The resulting vacancy-vacancy pair correlation function is shown in Figure 2.

The pair correlation function is also calculated using an analytic representation of the radial concentration profile using Eqs. (1) and (2). Coefficients in the analytic formula are determined from a multi-parameter least-square minimization of the difference between the pair correlation functions calculated from the defect maps and from the analytic expression. Good fits between the two pair correlation functions were obtained using

as the analytic representation of the defect concentration, where N is the number of defects in the cluster and λ is the mean radius. Figure 2 shows that the pair correlation functions calculated from the defect maps and from the analytic approximation (Eq. (5)) agree closely. Since Eq. (5) has a singularity at the origin, the defect concentration used in simulations was truncated at a value of 1 at. % to avoid unphysically large values at small radii and to represent rapid spontaneous recombination of defects when their separation approaches a lattice constant.

Equation (5) also gives good fits to pair correlation functions from defect maps for monoenergetic recoils. Figure 3 shows values of N and λ determined from such fits versus recoil energy. The curves in Figure 3 can be parameterized to determine values of N and λ as functions of recoil energy, which can be used to evaluate effective values of N and λ for irradiation conditions which produce a spectrum of recoil energies P(E) using

This simplification avoids the necessity of performing separate MARLOWE calculations and pair-correlation analysis fits for different neutron or recoil energy distributions.

Another significant approximation in the model is that the collision cascade is treated as a collection of point defects, whose properties, such as energies of formation and migration, are not influenced by neighboring defects. Thus, defect and carrier reactions in the central core region and in the lower concentration outer regions of the cluster are treated in the same way. In reality, the concentration of defects in the core is so high that this approximation may break down. Molecular dynamics simulations show that the lattice structure in the core region may be disordered or amorphous-like.⁵ Electronic properties of small dense defect clusters have been examined by density functional theory in GaAs⁶ and Si.⁷ These were found to have many electronic states and energy levels in the band gap and thus should contribute to carrier recombination. However, the contribution of the high-density core to overall carrier recombination can be limited by carrier transport into this region, in which case atomistic details are less influential. The significance of this approximation thus depends qualitatively on whether carrier capture or transport is the rate determining process.

MARLOWE and molecular-dynamics simulations show that the defects produced in the initial collision cascade include primarily interstitials, vacancies, and anti-site defects on the Ga and As sublattices, plus divacancies and di anti-site defects.⁸ Molecular-dynamics simulations of collision cascades were done using the LAMMPS code⁹ with bond-order-potentials for GaAs.¹⁰ The selection of defect species and reactions explicitly treated by the model is guided by calculations of defect properties in GaAs using DFT. These DFT calculations were done using the Socorro code¹¹ which uses a plane-wave basis set and the SeqQuest code^12–15 which uses a Gaussian basis set.

DFT determines the formation energies of various atomic configurations of the various charge states of defects. A defect in a given charge state may have multiple configurations with local energy minima with respect to small displacements, i.e., are metastable. The configuration with the lowest formation energy is referred to here as the stable configuration for that charge state. The stable configuration may be different for different charge states of a defect. The model explicitly tracks the concentrations of defects in their lowest-energy or stable state. Metastable configurations may influence carrier capture-emission reaction rates as discussed below.

The allowed charge states determine which carrier capture/emission reactions can occur, and their formation energies determine the relative populations of the various charge states in thermal equilibrium. Table I lists the formation energies for the lowest energy configuration, including lattice relaxation, of the various charge states of defect species included in the model. These values were obtained using the Socorro plane-wave DFT code with a newly developed bounds analysis approach.¹¹ These calculations included checks for cell-size convergence and Brillouin zone sampling. The methods used for these calculations are described in Ref. 11 for the case of anti-site defects. The same methods were used to obtain the values listed in Table I for the other defects. The model explicitly includes all primal defects: interstitials, vacancies, and anti-site defects of both types and divacancies. Carbon and silicon interstitials, produced by kick-out reactions between arsenic or gallium interstitials and dopants, may undergo further reactions and are therefore also explicitly included. The di anti-site defect was not included because DFT simulations show that no charge states other than the neutral state are stable and therefore this defect cannot participate in carrier recombination.^12,13 Table I includes a generic defect type intended to collectively represent immobile defect types not otherwise explicitly included, for example, products of defect reactions. The use of a generic defect type is motivated by computational simplification, to limit the number of species tracked and reactions included in the model. These generic defects were chosen to have +1, 0, and −1 charge states with formation energies to give energy levels at one third and two thirds of the energy gap for the 0/+1 and −1/0 levels, respectively. These were chosen to approximate the effect of charge states and energy levels for various defect species examined by DFT, which span the full range of the bandgap.

Defects may migrate by a thermally activated process, in which thermal fluctuations in their local atomic configuration lead to a higher energy saddle-point configuration from which the lattice relaxes to a minimum-energy defect configuration on an adjacent site. This leads to random-walk diffusion and to drift in electrostatic fields for charged defects. Migration barriers were determined from DFT calculations of formation energies of relaxed lattice configurations as the defect moves from one site to a neighboring equivalent site. Various pathways or intermediate configurations were explored to find the lowest migration energy. Of the primal defects initially produced in a collision cascade, only the arsenic and gallium interstitials are potentially mobile at room temperature. These may therefore diffuse and react with dopants and other defects to produce secondary reacted species. Kick-out reactions, in which an arsenic interstitial displaces a carbon p-type dopant from an arsenic lattice site, or a gallium interstitial displaces a silicon n-type dopant from a gallium site, produce carbon or silicon interstitials. DFT predicts that carbon and silicon interstitial atoms may also migrate and react to form tertiary products. The model therefore includes as potentially mobile species, arsenic, gallium, carbon, and silicon interstitials. All of the other defect species examined by DFT had large energy barriers for migration and are therefore not significantly mobile by thermally activated diffusion at room temperature. Table II lists the energy barriers for migration of the interstitial species obtained using the Socorro DFT code.¹¹ 

Thermally activated migration of defects can be described in terms of a diffusion coefficient

which is isotropic in a lattice with cubic symmetry. Values for D are used in evaluating defect reaction rates and the transport of mobile defects by drift and diffusion. Values for the diffusion prefactor D₀ are not accurately known, but can be estimated within about an order of magnitude from

for a cubic lattice, where d is the hop distance between neighboring sites, and ν is the attempt frequency. Using the optical phonon frequency as the attempt frequency, ν = 8×10¹² Hz (Ref. 16) and the spacing between like atoms for the hop distance $d=a0/2$ = 0.40 nm, Eq. (9) gives D₀ = 0.002 cm²/s. In view of uncertainties, we use D₀ = 10⁻³ cm²/s as a nominal value for the diffusion prefactor. Defects with migration energies below about 0.5 eV have sufficient mobility to migrate at a significant rate at room temperature. The migration energy is seen to vary significantly with charge state, giving rise to a strong dependence of the defect transport rate on the relative charge state populations. This leads in turn to acceleration of defect reaction rates under carrier injection, as will be discussed in Sec. IV.

DFT also shows that an additional athermal process for defect migration occurs for the arsenic interstitial.¹¹ This is the Bourgoin-Corbett mechanism^17,18 in which a change in charge state causes a change in the minimum-energy configuration such that the new configuration is the same as that of the initial charge state at the saddle point along the migration path. The interstitial can then migrate athermally by alternating capture of electrons and holes between the +1 and 0 charge state. The process is implemented approximately in the model as an additive term to the diffusion coefficients of As_I(+1) and As_I(0), using Eq. (9) with the frequency ν equal to the carrier capture rate per defect, and the hop distance d equal to the lattice constant times a scale factor with nominal value of 1. Thus, As_I(+1) makes one hop for each electron capture and As_I(0) makes one hop for each hole capture. This increased diffusion increases both transport and reaction rates for these two species.

Another possible mechanism for athermal or recombination-induced migration is the energy-kick process in which the electronic energy released by carrier capture converts to local atomic motion of sufficient magnitude and direction for the defect to cross a migration saddle point and move to a neighboring site.^19–21 This process is presently not implemented in the model because available information is not sufficient, though it could be included in a manner similar to that described above for the Bourgoin-Corbett mechanism.

Mobile interstitials react with immobile defects and dopants. Charge states and formation energies of the various reacting and product species, and hence the reaction enthalpies, were determined using the SeqQuest DFT code.^12–15 Table III lists the reactions between mobile interstitial species and the various primary immobile species examined by DFT. These DFT calculations identify the defect reaction network and which reactions are exothermic and hence likely to occur.

The model uses a standard drift-diffusion treatment for transport of both carriers and defects. The main influence of carrier transport in the present application occurs through drift in high fields where the carrier drift velocity saturates due to increased scattering by optical phonon and inter-valley transfer mechanisms. Under these conditions, the carrier mean-free path is small compared to the dimensions of a defect cluster justifying a drift-diffusion treatment of transport.

Continuity equations specify the rate of change of concentration n_i for the various species i being followed, which include conduction electrons, holes, dopants, and defects in their various charge states. For spherical geometry, these equations are of the form

where x is the radius. The reaction terms give the rate of removal or creation of species i through reactions with other species such as carrier capture/emission or defect reactions.

For transport by drift-diffusion, the flux of charge carriers can be expressed in terms of a gradient of the chemical potential ψ

where e is the elementary unit of charge and the mobility

is defined in terms of a drift velocity, v_drift in an electric field F.

For electrons and holes, the spatially dependent part of the chemical potential depends on the local electrostatic potential V and a concentration dependent term ψ_c

where the ψ_c are expressed relative to the corresponding band edge. The conduction electron chemical potential corresponds to the electron quasi-Fermi level within an additive constant. The hole chemical potential corresponds to the negative of the hole quasi-Fermi level within an additive constant.

When the concentration of electrons or holes is small compared to the density of states N_C for the conduction band and N_V for the valence band, Boltzmann statistics apply and ψ_c can be approximated as

where k is Boltzmann's constant and T is the temperature. In this case, Eq. (11) takes the familiar form

Here, we use a generalized form of the Einstein relation between diffusivity D and mobility μ for species with multiple charge q = Qe, where e is the elementary charge and Q is the charge state

When carrier concentrations are not small compared to the density of states, the carrier chemical potential is given by

where $F1/2−1$ is the inverse function of the Fermi-Dirac integral for a parabolic band

and Eqs. (14) and (15) no longer apply.

For mobile defects, Boltzmann statistics is applied. This results in an analog to Eq. (15) for the flux where mobility is replaced by the diffusivity using Eq. (16) including the charge q of the defect

The electric potential V is determined from the charge density ρ by solving Poisson's equation, which for spherical geometry is

where ε is the dielectric constant of the material.

At high field, the carrier drift velocity saturates. The dependence of carrier mobility on electric field strength |F| is approximated by

for holes, where μ_o is the mobility at low field and v_sat is the saturation velocity at high field.²² For electrons in GaAs, the field dependence of mobility is calculated from

to include the transferred electron effect, where F₀ is a saturation field strength.²³ The low field mobilities are input parameters whose value depends on doping concentration and temperature.

Evaluation of the rates of carrier capture and emission by defects is an important aspect of the model, since these rates influence the net carrier recombination rate. Moreover, the carrier reaction rates determine the population of defect charge states, and since the activation energy for defect migration depends on their charge state, their diffusion, and reaction rates are also affected.

This section discusses the formulation of rate equations for carrier capture and emission reactions. We first formulate the capture and emission of carriers by defects in terms of one-step conversions between charge states in their lowest-energy configurations listed in Table I, neglecting any higher-energy intermediate configurations. Reactions for some particularly influential defects, where there is major reconfiguration and more detailed information on the reaction path has been obtained from DFT, are formulated as a two-step process. In both cases, the rate equations are formulated in terms of a cross section for carrier capture. Values for these cross sections are discussed. We conclude the section with a discussion of the enhancement of carrier reaction rates by band-to-trap tunneling in the electrostatic field produced by charging of defects.

The reaction terms in Eq. (10) remove or create carriers and defect species. Reactions involving charge carriers are treated as reversible and are formulated to give correct thermodynamic equilibrium in detailed balance. The electron-hole reactions included are direct recombination and Auger recombination

where A is the coefficient for direct recombination and B_e and B_h are coefficients for Auger recombination. The factor $(nenh)0$ is the product of electron and hole concentrations in thermal equilibrium under the constraint of constant local charge density.

The model includes rate equations for the capture and emission of charge carriers by allowed charge states of the defects listed in Table I giving 56 carrier-defect reactions. In the case of electrons, such reactions are of the form

In the forward or exothermic direction, a conduction electron is captured by a defect of type Z in charge state Q, while in the reverse reaction an electron is emitted from a defect of type Z in charge state Q − 1. ΔE^e(Z_Q−1) is the positive difference between the formation energies of the initial and final states

where E_f0(Z_Q) and E_f0(Z_Q−1) are the formation energies, referenced to the valence band edge, of the defect in the two charge states in their lowest energy configuration, and E_g is the band gap energy.

For reactions in which the symmetry of the defect and local bonding configuration is maintained and the reconfiguration is modest, akin to a breathing mode relaxation, we use a single-step process in which the capture rate is proportional to a cross section for carrier capture σ times a carrier flux equal to the product of the carrier concentration and thermal velocity. The corresponding carrier emission reaction is formulated to give the correct concentrations of the two charge states in thermal equilibrium. The rate equation for the above reaction, including both capture and emission reactions, is

for Fermi-Dirac carrier statistics. We use the notation U{X,Y} to denote the rate of reaction between species X and Y. n_e, [Z_Q], and [Z_Q−1] are the concentrations of the three species involved, v_th^e is the electron thermal velocity, $σZQe$ is the cross section for capture of an electron by defect Z_Q, and $ψce$ is the electron chemical potential defined in Eq. (17).

The corresponding reaction equations for hole capture and emission are

where ΔE^h(Z_Q+1) is the energy to emit a hole from defect Z_Q+1 to the valence band, and

Note that ΔE^h(Z_Q+1) + ΔE^e(Z_Q) = E_g.

DFT shows that some charge-change reactions occur between states with very different structures, with different symmetry or type of neighboring atoms. Figure 4 shows the formation energies relative to the valence band edge (E_VB) from DFT calculations, for the configurations and charge states of the arsenic, gallium and carbon interstitials, which are most consequential for the present model. Vertical arrows indicate a charge change reaction by emission (or capture) of a hole to (or from) the valence band between states with similar configurations. Diagonal arrows indicate reconfiguration between states of the same charge. Corresponding diagrams for electron capture and emission can be obtained by plotting formation energies E_CB relative to the conduction band edge instead of the valence band edge, where E_CB = E_VB + QE_g.

Some charge-change reactions occur via metastable intermediate states for which there is an energy barrier for reconfiguration. Examples of this are shown in Figure 4(b) for the Ga_I{T_d:a, + 2} ⇔ Ga_I{T_d:g, + 1} + h reaction which has a reconfiguration barrier of 0.93eV, and in Figure 4(c) for the C_I{C_2V:001g, + 1} ⇔ C_I{C_1h:pAs001a, 0} + h reaction which has a reconfiguration barrier of 0.35 eV. The above notation of defect state denotes: defect type{configuration, charge state}.

Large relaxation energies and energy barriers for reconfiguration have strong effects on carrier capture and emission rates which are not adequately captured by a single-step process between initial and final states as described in Eqs. (26) and (28). For these cases, we treat the reaction as a two-step process of change in charge between states of similar configurations and reconfiguration between states of the same charge. The rate equations are formulated without explicitly following the population of additional intermediate states by considering a steady-state condition. These rate equations use energies for reconfiguration and carrier emission, and the energy barriers for reconfiguration calculated by DFT, which are listed in Table IV and illustrated in Figures 4(a)–4(c).

Equations for reaction rates for the above two-step process are derived as follows. Let index i = 1 represent the initial, stable state of the defect that will emit a conduction electron or hole; i = 3, the final stable state that will ensue; and i = 2, the intermediate state. The atomic reconfiguration may precede or follow the emission. The reverse process involves carrier capture. In terms of concentrations C_i (units of inverse volume) and i-to-j reaction-rate coefficients R_ij (inverse time), one has

Consider the approach to a steady-state reaction condition while the stable-state concentrations C₁ and C₃ are fixed. Rewriting Eq. (29) gives

which exhibits the time constant and the asymptotic value of C₂. The steady-state reaction rate per unit volume is then

Equation (31) has two limitations for the description of system behavior; transient phenomena on the scale of the time constant of Eq. (30) are not fully captured; and, the concentrations of defects in the intermediate states do not explicitly appear in the model.

The above equations have been implemented as follows for the two circumstances that arise during reversible defect charge exchanges accompanied by atomic reconfiguration:

1. endothermic reconfiguration of the defect at constant charge followed by endothermic carrier emission at fixed configuration, accompanied by the two corresponding reverse reactions, denoted by [RE] for reconfiguration-emission;

2. endothermic carrier emission followed by exothermic reconfiguration, and the reverse reactions, denoted by [ER] for emission-reconfiguration.

Exothermic reconfiguration at constant charge with a positive activation barrier ΔE_bar (expressed relative to the starting configuration) is assumed to occur with a prefactor ν₀.

One then has

Detailed balance in thermodynamic equilibrium then gives for the reverse, endothermic reconfiguration

where ΔE_cfg is the positive change in configuration energy. The prefactor ν₀ is a characteristic frequency for lattice vibration which we approximate as the optical phonon frequency.¹⁶ 

The rate coefficient for exothermic carrier capture is

where n is the concentration of carriers (electrons or holes). Detailed balance then gives for the reverse, emission event

for Fermi-Dirac carrier statistics, where ΔE_em is the positive carrier emission energy and ψ is the carrier chemical potential given by Eq. (17). Table IV lists the reactions for which the two-step process is used in the model, and the values for the reaction energies: ΔE = ΔE_em + ΔE_cfg is the total difference in formation energy, ΔE_em is the activation energy for carrier emission at constant configuration, ΔE_cfg is the energy of reconfiguration at constant charge state, and ΔE_bar is the energy barrier for reconfiguration.

For some charge-change reactions, there is no metastable intermediate state, for example, as shown in Figure 4(a) for the arsenic interstitial. DFT shows that bound states exist for the +1 charge state in all three configurations, but the As_I{T_d:g, + 1} and As_I{C_2v:110a,+1} configurations are both dynamically unstable and relax with no barrier to the lower energy As_I{C_1h:p001g,+1} configuration. For reactions of this kind, we apply the two-step formulation using an intermediate state which is lattice-relaxed under the constraint that the symmetry is the same for both charge states. This allows breathing-mode type relaxations to occur, similar to those in the case of one-step reactions or two-step reactions with a metastable intermediate state. However, we recognize the possibility that such reactions might also occur via other configurations with lower formation energy. We therefore consider this approximation to give an upper bound on the energy barrier and a lower bound on the rate for such reactions. The single-step process with no intermediate state would give an upper bound on the reaction rate.

Figure 4(a) also illustrates the Bourgoin-Corbett process for athermal diffusion, which occurs by alternating capture of electrons and holes:

• As_I{C_1h:p001g, +1}+e → As_I{C_1h:p001g, 0} relaxes to As_I{C_2v:110a,0}, and

• As_I{C_2v:110a, 0}+h → As_I{C_2v:110a, +1} relaxes to As_I{C_1h:p001g, +1}.

Equations (26)–(35) for the rates of carrier capture and emission are formulated in terms of a cross section σ for carrier capture between states with similar configurations. Values of capture cross sections for various defects in GaAs have been reported from experiments.^24–28 This information comes mainly from measurements of carrier emission rates by deep-level transient spectroscopy (DLTS) in which the identity and structure of the defect is unknown in most cases. To explore behavior of our model, we have provisionally used nominal values for carrier capture cross sections based on these experimental values. Nominal values for capture cross sections used here are 10⁻¹⁶ cm² for capture by a neutral defect, 10⁻¹⁵ cm² for capture by a charged defect where the Coulomb interaction is attractive, and 10⁻¹⁷ cm² for capture by a charged defect where the Coulomb interaction is repulsive. For a two-step reaction process, the rate will be determined by the slower of the two processes. For reaction by carrier capture followed by exothermic reconfiguration with no barrier, the carrier capture step will generally be slower than the reconfiguration. In this case, the two-step and one-step processes give the same reaction rate.

We anticipate that more accurate values for carrier capture cross sections will become available. Models for carrier capture have been developed based on a multi-phonon process,^26–32 which are often formulated in terms of a configuration coordinate and parabolic potentials. These models are appropriate for charge-change reactions where the lattice relaxation is of the breathing-mode type. The relaxation couples the electronic energy to phonons. The carrier capture rate becomes temperature dependent and can be described in terms of a capture cross section with an activation energy which depends on the binding energy of the carrier to the defect, lattice relaxation energy, and phonon energy. DLTS measurements provide values for these parameters and for capture cross sections and their dependence on temperature,²⁸ however, the defect types corresponding to the DLTS peaks are not known in most cases. It may be possible to associate DLTS data with specific defect types using defect properties from DFT calculations such as carrier emission energies, activation energy for thermal migration, and lattice relaxation energies. Furthermore, DFT-based methods for calculation of carrier capture rates from first principles, recently developed for GaN,^33,34 may also be applied to defects in GaAs.

High electric fields may greatly increase the rate of carrier capture and emission by traps through band-to-trap tunneling with energy coupling to phonons.^28–31 Such high fields arise in depletion regions near junctions in devices. Charging of clustered defects can also produce large electric fields, particularly in heavily doped neutral regions as typically used in the base of HBTs. In this band-to-trap tunneling process, field-induced band-bending lowers the relative energy of band states in the vicinity of a trap, both increasing the carrier population of these states and reducing the energy that must be exchanged with the phonons during the tunneling-enabled capture. To enable exploration of its potential impact on carrier recombination at clustered damage, we have implemented band-to-trap tunneling in our simulation, based on a formulation for a multi-phonon carrier capture process described by Schenk.³⁰ Two extensions from Schenk's formulation are necessary due to the fact that field-enhanced recombination rates combined with very high defect concentrations and transport-limited flow of carriers into the cluster can reduce carrier concentrations within the cluster by orders of magnitude from equilibrium values. In our formulation, the relevant carrier concentrations are therefore extracted as a function of location from the model simulation over the range of tunneling instead of assuming that the carrier concentrations are in equilibrium; and carrier removal from the bands is appropriately delocalized rather than being taken to occur entirely at the position of the trap. Also, detailed balance is applied to determine the rate of carrier emission from the traps for correct equilibrium. The influence of band-to-trap tunneling on carrier recombination will be discussed in Sec. IV. Details of the formulation and implementation of band-to-trap tunneling will be given in a separate publication.

The sink term for defect-defect and defect-dopant reactions of the type X + Y → Z, between species of type X and Y producing a product of type Z is evaluated using

One of the reacting species must be mobile, i.e., an interstitial. The defect reaction rate is proportional to the product of the concentrations [X] and [Y] of the two reactants, the diffusion coefficient D of the mobile defect, and an effective reaction radius R. In some cases, there is no product defect, for example, in the annihilation reaction between an interstitial and vacancy of corresponding type. Defect reactions are configured to conserve total charge. When a product defect is produced, the charge state of the product is the sum of the charge states of the reactants, provided that there is a stable configuration of the product defect as determined by DFT. When the reaction would produce a product defect in a charge state that is not stable, the reaction is taken to produce the defect in the nearest stable charge state, plus electrons or holes to conserve charge. Examination of reaction enthalpies, based on defect formation energies determined by DFT, shows that most defect reactions are sufficiently exothermic to be essentially irreversible at room temperature. Therefore, inverse dissociation reactions are presently not included in the model (Eq. (36)), although such reactions could be added if needed. The number of defect-defect and defect-dopant reactions included in the model is 286.

The reaction radius R in Eq. (36) is evaluated as follows. If either reacting species is neutral, the reaction radius is equal to the lattice constant a₀ = 0.565 nm for GaAs. When both species are charged with the same sign the reaction is not included, under the assumption that Coulomb repulsion inhibits the reaction. When the two species are charged with opposite sign, the reaction radius is the distance at which the energy of the screened electrostatic interaction equals the thermal energy

with Debye screening length $λD=εkT/e2n$ from mobile charge with concentration n and where Q₁ and Q₂ are the charge states of the reactants.

The code uses an analytic approximation for the root of Eq. (37) to evaluate the reaction radius

where |Q₁ Q₂| R₀ is the reaction radius without screening, and R₀ = e²/(4 π ε k T) = 4.25 nm for GaAs at 300 K. Where the defect concentration is high, the reaction radius calculated from Eq. (38) may be greater than the mean distance between defects, approximated as R_def = N_def^−1/3, where N_def is the defect concentration. The reaction radius is therefore reduced so that it does not exceed half the distance between defects

Numerical solution of Eq. (10) is accomplished through spatial discretization on a radial mesh x_j which may have variable step size. This converts the partial differential equations to a coupled system of ordinary differential equations (ODEs), which are then integrated forward in time from a specified initial condition. The extreme numerical stiffness of the problem and the highly sensitive interdependency of variables require a suitably advanced solver. We employed the ODE solver DDEBDF.³⁵ Solutions are formulated to produce correct thermodynamic equilibrium and to conserve the number of particles and charge.

Discretization of Eq. (10) gives

for the time derivative of concentrations of species i at radial mesh index j > 1, where $Φji$ is the flux from radius x_j−1 to x_j. Here, we use the notation that the superscript denotes the species index and the subscript denotes the radial mesh index. At the origin, j = 1 and x₁ = 0 and

Numerical evaluation of the flux is accomplished by a generalization of the approximation originally described by Scharfetter and Gummel³⁶ for Boltzmann carrier statistics

where

V_j is the electrostatic potential at the mesh points and q = Qe is the charge (with charge state Q = −1 for electrons, +1 for holes). Equation (42) is exact for steady-state flow under uniform field, yields thermodynamic equilibrium correctly at zero flow, and approaches the correct derivative form as Δx goes to zero. The flux of mobile defects with charge q ≠ 0 is calculated using Eq. (42) with an equivalent mobility derived from the diffusivity using Eq. (16). For neutral defects, the flux is evaluated from

where $Δnji=nji−nj−1i$⁠.

Equation (42) is used to evaluate the fluxes of charged defects and of carriers when their concentrations are small compared to their effective density of states. When this condition is not satisfied, Fermi-Dirac statistics must be used for the carriers. In this case, the carrier flux is calculated using a modified form of Eq. (42)

where Q = −1 for electrons and Q = +1 for holes and

are evaluated from the carrier concentrations using Eq. (17). The form of Eq. (45) was chosen to satisfy three criteria: reduction to Eq. (42) when carrier concentrations are small compared to the density of states; correct thermodynamic equilibrium at zero flow; and reduction to the correct derivative form as Δx goes to zero. Equation (45) is used only at the higher concentrations where it is required; Eq. (42) leads to greater stability in the numerical solution at small concentrations, where zero and even negative concentration values can arise temporarily as the calculation proceeds. Reaction terms in Eqs. (40) and (41) are calculated by direct application at the mesh points x_j of Equations: (23) for direct and Auger recombination, (26) and (28) for carrier capture/emission by defects for the one-step process, and (31) for the two-step process, and (36) for defect-defect reactions. Although Eqs. (27) and (28) are mathematically correct for all carrier concentrations >0, computational considerations make it advantageous to evaluate the reaction rates using corresponding equations for Boltzmann carrier statistics

for electrons, and

for holes, when the carrier concentration is small relative to the effective density of states, and using Eqs. (27) and (28) for Fermi-Dirac carrier statistics otherwise.

Poisson's equation is solved at each time step. The numerical method used is to advance the spatial integration of Eq. (20) from the center of the sphere to the outer boundary using the exact solution for the case where the charge density varies linearly between mesh points, ρ/ε = γ + βρ, with the boundary conditions that the field and potential are zero at the center. This integration uses the same spatial discretization as that used for solving the transport equations. The discretized equations for this case are

At the center (i = 1) x₁, V₁, F₁, a₁, and b₁ all have values of zero.

These equations give the potential V_i and field F_i at the mesh points x_i. Evaluation of the fluxes (Eqs. (42) and (45)) uses these potentials at the mesh points, but employs carrier mobilities corresponding to the field midway between mesh points, approximated as

The Poisson solver allows two options for charge compensation:

• no compensation,

• subtract uniform background charge to give zero combined charge density at the boundary.

The latter option is used for simulations of recombination in a depletion region where charge on ionized dopants would otherwise give a spurious radial electric field not present in a device.

The code calculates carrier chemical potentials from their concentrations according to Eq. (17), using the Joyce-Dixon approximation³⁷ for the inverse Fermi integral for parabolic bands when x < 7.5, and an expression suggested by Blakemore³⁸ 

when x > 8.5, where x = n^e/N_C for electrons and x = n^h/N_V for holes. The gap between 7.5 < x < 8.5 is bridged using a pair of connected parabolas chosen to maintain a continuous first derivative.

At the outer boundary of the sphere, a reflecting boundary condition (flux = 0) is always used for mobile defects. The code provides three options for the boundary condition on carrier concentrations at the outer boundary

• concentration of electrons and holes are set to the thermodynamic equilibrium value,

• no flow of carriers through the boundary, i.e., reflecting boundary and net contained charge is conserved,

• concentrations of electrons and holes are set to constant specified values.

Because we are using an ODE solver, values for the time derivative of concentrations at the boundary must be specified rather than values of concentrations. Thus, boundary conditions (a) and (c) are implemented by setting the time derivatives at the boundary to

where $τ=Δx2/D$ is a time constant for diffusion between mesh points with diffusion coefficient $D=μkT/e$⁠. This is a source/sink term for carriers at the boundary which causes the value of the carrier concentration at the boundary to rapidly approach the specified value n₀, without destabilizing the ODE solver.

The boundary conditions on the solution of Poisson's equation are that the electrostatic potential and electric field are both zero at zero radius.

The initial condition used in the model is that the concentrations of electrons and holes are equal to the concentrations of n and p-type dopants, respectively, the dopants are ionized, and the defects are initially in their neutral charge state. The initial radial distribution of defects within a cluster is obtained from a pair correlation function analysis of defect maps calculated using the binary collision code MARLOWE as described in Sec. II A. For the simulations discussed below, the MARLOWE calculations used a distribution of recoil energies for neutrons with an energy spectrum for a fast-burst reactor,³ as illustrate in Figure 2. The defect concentration versus radius is calculated using Eq. (5) truncated at a concentration of 1 atomic percent or 4.4 × 10²⁰/cm³, with parameter values of N = 550 for the number of defects per cluster and λ = 11.6 nm for the mean radius for this irradiation condition.

Although concentrations of vacancies in the central core region of a cluster is typically slightly higher than that of interstitials, in the simulations discussed below the same radial concentration profiles were used for vacancies and interstitials since the concentrations used in this core region are truncated and therefore approximate anyway. Also, the same radial profiles were used for gallium and arsenic vacancies and interstitials. The initial concentrations of anti-site defects and divacancies were taken to be zero in the simulations discussed below. More realistic values for the concentration of antisite defects and divacancies could be obtained from molecular dynamics and MARLOWE simulations but these were not included in simulations discussed here. The effect of their inclusion should mainly be to increase carrier capture rates, qualitatively similar to that of vacancies.

Properties of GaAs used in the model apart from radiation effects are widely documented^16,39–41 and are not discussed in detail here.

This section presents results from simulations for conditions relevant to operation of an Npn HBT and chosen to illustrate important features of the model and their implications for radiation effects. First, we consider the p-type base since carrier recombination in this region will reduce the forward-active gain. Figure 5 shows a plot of the total number of the various types of defects tracked by the model as a function of time from their creation. The time regime of interest is microseconds to seconds. The initial defect types are Ga and As interstitials and vacancies in equal numbers with a radial distribution for clustered damage produced by neutrons from a fast-burst reactor. For this simulation the concentration of carbon dopants was 5 × 10¹⁹/cm³ and the carrier concentrations at the boundary were held constant at 5 × 10¹⁹/cm³ for holes, corresponding to the neutral or field-free region of the base, and 1 × 10¹⁵/cm³ for electrons corresponding to the minority carrier concentration in the base of an Npn HBT with high operating current density as determined from device simulations. For this doping, values of 2240 cm²/V s for electrons and 30 cm²/V s for holes were used for the low-field mobilities and 7.7 × 10⁶ cm/s for the saturation velocity for both electrons and holes.^16,39–41 The maximum radius, i.e., radius of the outer boundary of the simulation was chosen to be 100 nm. At this value, the defect concentrations and electric field are very small, and the carrier flow into the cluster is insensitive to the boundary radius value.

Around 10⁻⁷ s As interstitials react, mainly with carbon dopants to produce carbon interstitials. A few As interstitials also annihilate As vacancies and react with Ga vacancies to produce As(Ga) antisite defects. Around 10⁻⁴ s the carbon interstitials react, mainly with other carbon dopants producing an immobile reaction product. Gallium interstitials do not react significantly. Figure 5 also shows the equal flows of electrons and holes into the cluster, indicating the combined rate of carrier recombination. This simulation used nominal values for cross sections for carrier capture as discussed in Sec. II D 1.

Figure 6 shows results from a simulation similar to that of Figure 5, with the same number of defects now distributed uniformly instead of clustered. Uniform damage is produced by irradiation with MeV electrons for which recoil energies are smaller and create single Frenkel pairs instead of displacement cascades. The defect reactions are similar for clustered and uniform damage, the main difference being the absence of reactions with vacancies for uniform damage. This difference is due to the fact that with uniform damage the defect concentrations are much lower than the concentration of carbon dopants, whereas with clustered damage the defect concentrations exceed the carbon concentration in the core of the cluster. Another difference is that the carrier recombination rate is about seven times larger for clustered damage than for uniform damage. This difference is due to the influence of band-bending on the concentration of carriers in the core of the cluster and on the flow of carriers into the cluster. The large hole concentration causes the defects to acquire a net positive charge. This localized charge changes the local potential and produces a high electric field in the cluster as shown in Figure 7. This increases the concentration of electrons in the cluster core and increases the drift term in their transport into the cluster. Figure 7 also shows the carrier flow versus radius. Carrier recombination occurs predominantly at radii less than 10 nm.

As discussed in Sec. II D 1, the values for carrier capture cross sections are poorly know. Therefore, we have examined the sensitivity of the model to variations in values of carrier capture cross sections. Increasing all carrier capture cross sections by a factor of 10 increased the carrier flow by a factor of 10 for uniform damage, as one might expect, but for clustered damage the carrier flow increased by only a factor of 1.6. This indicates that for clustered damage with nominal values for carrier capture cross sections, the carrier flow, i.e., recombination, is limited primarily by transport rather than by the rate of capture in heavily doped neutral material. This conclusion was confirmed by simulations in which the carrier mobilities were changed.

Carrier reaction rates can also influence defect reaction rates, through their effect on the population of the various defect charge states, due to the fact that defect mobility depends on charge state. A striking example of this is the carbon interstitial. In the simulation of Figure 5, elimination of the 0.35 eV barrier to reconfiguration in the (+1/0) carrier capture/emission reaction (illustrated in Figure 4(c)) caused the less mobile +2 charge state to become the most populated state, instead of the more mobile neutral charge state, greatly slowing the reaction rate for the carbon interstitial.

Figure 7 shows that at radii where carrier recombination is high the electric field is also high. Previous studies using DLTS show that carrier capture and emission by defects in GaAs are greatly enhanced in fields above 0.1 MV/cm by the effect of the field on band-to-trap tunneling.²⁸ Consequently, we have implemented this process into the cluster model and examined its influence on carrier recombination. Simulations were done for a cluster of immobile defects in heavily carbon doped material with band-to-trap tunneling properties of the strongest of the traps found in the earlier DLTS study.²⁸ The main effects of adding the tunneling were: (i) the carrier flow into the cluster increased 3%, the radius where carriers recombine increased to 10 to 15 nm, and (ii) the concentrations of electrons and holes decreased by several orders of magnitude at radii less than 8 nm. These effects are all consistent with the above conclusion that carrier recombination at clustered defects in the neutral base is flow-limited. Apart from these simulations done specifically to examine effects of band-to-trap tunneling, this process was not included in other simulations discussed here.

Figure 8 shows the radial concentration profile of defects and carriers at 1 s for the simulation of Figure 5 with clustered defects. Also shown are the initial concentration profiles of As and Ga interstitials and vacancies used in the simulation. The main defects remaining at late times are Ga interstitials, Ga and As vacancies, As(Ga) antisites in the core and the generic immobile reaction product (dmg) in the wing produced by carbon interstitials reacting with carbon dopants.

Figure 9 shows the number of defects versus time from a simulation similar to that of Figure 5 with clustered damage, but for equilibrium, corresponding to conditions in the base of a HBT with no bias applied and no current flowing. The reaction of As interstitials is about 2000 times slower due to absence of recombination-induced migration by the Bourgoin Corbett process. The reaction of carbon interstitials is also much slower because the population of the mobile neutral charge state is reduced and the less mobile C_I(+2) is now the most populated charge state. This comparison demonstrates the feature of the model that forward active operation of a HBT greatly accelerates the rate of defect reactions. This is consistent with experimental evidence that gain degradation from displacement damage recovers faster when current is flowing through the device,⁴² a process sometimes referred to as current-induced annealing.

In silicon transistors, carrier recombination in the depletion region of the emitter-base junction often contributes significantly to base current and gain degradation from defects. However, in HBTs, the emitter has a wider bandgap than the base. Device simulations show that this wider bandgap depresses the minority carrier concentration in the depletion region, with the result that when defect concentrations and carrier capture cross sections are uniform throughout the device, recombination in the depletion region is small compared to recombination in the neutral base. Figure 10 shows the number of defects and carrier flow versus time from a simulation for conditions chosen to represent the depletion region of a HBT. In this simulation the carbon dopant concentration was 5 × 10¹⁷/cm², the electron and hole concentrations were 10¹² and 5 × 10¹³ /cm³ respectively. The carrier flow into the cluster is several orders of magnitude lower than for a defect cluster in the neutral base shown in Figure 5. Because the carrier flow is low, it is no longer transport limited and therefore depends more directly on the carrier capture cross sections. Furthermore, due to the lower carrier densities, the defect charging, band bending and associated electric field are much lower than for a defect cluster in the neutral base and the region of recombination extends to larger radii. Thus, the influence of defect clustering on carrier recombination is smaller in the depletion region than in the neutral base. However, the device field in a depletion region of a junction can be large enough to increase carrier recombination through field-enhanced band-to-trap tunneling. Although we have implemented field-enhanced band-to-trap tunneling in our model for carrier recombination in a defect cluster in bulk material, the situation is different in the depletion region of a HBT where the close proximity of the heterojunction is likely to influence band-to-trap tunneling. Simulations of reduction in gain of a HBT by defects require implementation of field-enhanced recombination by band-to-trap tunneling in a device model, including the influence of the device field in depletion regions.

Simulations were also done for conditions representing the neutral base of a Pnp HBT; with a concentration of silicon n-type dopant of 5 × 10¹⁸/cm² and carrier concentrations at the boundary of 5 × 10¹⁸/cm² for electrons and 10¹⁵/cm³ for holes. In this case the As interstitials react between 10⁻⁶ and 10⁻² s, their dominant charge states are 0 and −1, and their rate of reaction is slower with lower minority carrier concentration. The Ga interstitials do not react so formation of silicon interstitials by the kick-out reaction Ga_I + Si(Ga)-> Si_I, is insignificant.

The main objective of this work was to develop a model to simulate carrier recombination in clustered damage produced by atomic collisions in GaAs, thereby providing a basis for predicting the transient response of HBTs to pulsed neutron damage. The model incorporates the physics of charge carriers and defects in GaAs and their reactions, and the influence of defect clustering. Information on properties of defects was obtained mainly from recent DFT simulations.

Results from simulations for various conditions were discussed, focusing mainly on heavily carbon doped p-type GaAs corresponding to the neutral base of an Npn HBT. In heavily doped neutral material, clustering of the defects has a strong effect on carrier recombination. Charging of the defects changes the local electric potential creating a large electric field within the cluster. This increases the concentration of minority carriers within the cluster and increases their transport into the cluster by field drift, both of which increase the rate of carrier recombination. On the other hand, with nominal values for carrier capture cross sections, the carrier flow into the cluster, i.e., the recombination rate, was found to be limited by transport of minority carriers rather than by the rate of capture in heavily doped neutral material. In contrast, for conditions representing the depletion region at the emitter-base junction of a HBT, the influence of clustering on carrier recombination is weaker. In this case, the carrier concentrations are lower and the localized defect charge, and hence the cluster electric field, are smaller. Consequently, rates of carrier recombination at defects should be much lower in the depletion region of a junction than in the neutral base. However, high device fields in depletion regions might increase carrier recombination at defects by enhancement of band-to-trap tunneling, and the impact of this on device operation remains to be explored.

Another important aspect of the model is the time dependence of carrier recombination due to the migration and reaction of defects. Simulations for conditions in the heavily carbon-doped base of an Npn HBT show that arsenic interstitials react with carbon dopants very quickly producing carbon interstitials by a kick-out reaction. These carbon interstitials also migrate and react more slowly, mainly with carbon dopants to form an immobile defect complex. In this case, the defect reactions are not greatly different for defects which are clustered versus uniformly distributed, possibly due to the fact that the defect concentration is below the carbon concentration except in the inner core of the cluster. In simulations, higher concentrations of minority carriers greatly increase rates of these defect reactions, consistent with experimental evidence that gain degradation of a HBT from displacement damage recovers more quickly with current flowing through the device. Mechanisms for recombination-driven annealing were identified as (i) Bourgoin Corbett process for As interstitials and (ii) an increase in population of more mobile charge states for carbon interstitials.

Capture-emission reactions between carriers and defects are described in terms of carrier capture cross sections, whose values affect both carrier recombination and defect reaction rates. Carrier capture rates are strongly influenced by reconfiguration of the defect structure associated with changes in charge state. DFT has provided new information on this effect by determining formation energies for multiple configurations of the various charge states for defects in GaAs. This information enabled a first approximation of the influence of defect reconfiguration on carrier capture rate to be incorporated in the model as a two-step process for a few reactions where the reconfiguration is large. Previous work indicates that more generally, lattice relaxation gives rise to a temperature dependence to the capture cross section, which is often approximated as an activation energy for carrier capture.^26–32 Recently developed DFT-based computational methods³³ are being applied to obtain more quantitative information on this temperature dependence which should reduce uncertainty in carrier reaction rates in the model.

Comparisons between the model and experiments are only qualitative at present. For example, the model reproduces the current-induced recovery of gain observed in irradiated HBTs.⁴² Quantitative comparisons are more easily made with data from experiments on bulk material, such as transient carrier lifetime after pulsed irradiation used previously for silicon.¹ Such data are not available for GaAs. Experimental information is available on the transient recovery of gain of irradiated HBTs, however, because of the extreme variation in conditions within the device, quantitative comparison to such data requires integration of the carrier recombination model with a device model. Such an integrated model will realize the goal to provide a capability for predicting the transient response of HBTs to displacement damage.

The authors thank Alan Wright for his DFT calculations of defect formation energies and of pathways and energy barriers for migration and charge-change reactions, and for his permission to use results from those calculations prior to their publication. The authors thank Peter Schultz for his DFT calculations of energies for defect formation and reactions. The authors acknowledge valuable contributions from Philip Cooper and Brian Hehr in providing defect maps for neutron irradiated GaAs and from Stephen Foiles for Molecular Dynamics simulations of collision cascades in GaAs, which together were used to determine a representative radial defect concentration profile for neutron damage. The authors acknowledge valuable discussions with David Lang, Andreas Schenk, and Normand Modine on the topic of carrier capture by defects. This work was conducted at Sandia National Laboratories under the auspices of the U.S. Department of Energy's National Nuclear Security Administration.

Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under Contract No. DE-AC04-94AL85000.