Non Local Impact Ionization Effects in Semiconductor Devices

Impact ionization processes define the breakdown characteristics of semiconductor devices. 
An accurate description of such phenomenon, however, is limited to very sophisticated device 
simulators such as Monte Carlo. A new physical model for the impact ionization process is 
presented, which accounts for dead space effects and high energy carrier transport at a Drift 
Diffusion level. Such model allows to define universal impact ionization coefficients which 
are device-geometry independent. By using available experimental data these parameters 
have been calculated for In0.53Ga0.47As.


INTRODUCTION
Device simulation has become a fundamental part in the development and optimization of electronic devices.Next to standard approaches such as the Drift Diffusion method ], which has lead to several com- mercial software packages, newer techniques are also available, ranging from hydrodynamic approaches to particle-based simulations, such as the Monte Carlo or the Cellular Automaton methods [2,3].
In general, one pays the enhanced physical content of the latter approaches with a heavier computational load.Thus is particularly true in the presence of high electric fields.While physical simulators require an accurate description of high energy effects, Drift Dif- fusion methods can rely on simplified phenomenological approaches.
The aim of this paper is to show the inadequacy of some of these approaches, by focusing on the process which is the main responsible for the breakdown of devices, namely impact ionization (II).Furthermore, we will indicate improvements which allow us to keep the speed advantages of Drift Diffusion methods while incorporating the main physical features of the II process.

IMPACT IONIZATION IN DRIFT DIFFUSION AND MONTE CARLO SIMULATORS
II is the phenomenon by which a hot carrier colliding with a valence electron creates a new pair of carriers, both available for conduction.A careful model of II must account for the high energy behaviour of carri- ers, since the threshold for the process is at least equal to the energy gap of the semiconductor under consid- eration.The best treatment in terms of physical accu- racy is provided by the Monte Carlo method, which Corresponding author.Tel: +39 6 72594469.Fax: +39 6 2020519.E-mail: lugli@roma2.infn.it291 typically includes II by associating to each carrier ionization probability per unit time, dependent on the carrier status (energy, velocity, momentum etc.).Sev- eral degrees of depth are possible, based on simplified bands and phenomenological II rates, or on complete bands and microscopical II rates.The latter approach has been applied mainly to bulk semiconductors, due to its numerical complexity, while the former one can be applied to devices.Recently, a Monte Carlo analy- sis has been presented of an AlGaAs/GaAs Heterojunction Bipolar Transistor (HBT) [4,5], which perfectly reproduced available experimental results both on the multiplication factor as well as on the electroluminescence spectra of the device in a near- breakdown regime [6,7].There, a three valley non- parabolic model was implemented for electrons and holes, together with the Kane model for II.The main conclusions of that work have been 1.The average energy and the ionization coefficients reach their maximum not at the base-collector junction, where the electric field reaches its maxi- mum, but rather inside the collector.Such effect is referred to as "dead space" effect [8][9][10][11].

Electrons and holes gain considerable energy in
the collector due to the presence of very high elec- tric fields.As a result, the electron and hole distri- bution finctions are very hot, leading to strong ionization processes and to radiative transitions within the conduction and the valence band, respectively [6], responsible for the observed elec- troluminescence.
Although the Monte Carlo simulation is physically very accurate, it is also extremely time consuming.It would therefore be preferable, in the context of device modeling (and especially in the presence of ionization phenonena), to use faster numerical tools such as those based on the Drift Diffusion algorithm.There, the following equations are used dp.
Jp qplJpF qDp--x dn Jn qnlJrtF + qDn--r- ax (2) (3) dx q(G-R) dx In Eqs.(3), II is accounted for by means of the gener- ation term, G.Unfortunately such approaches do not correctly describe hot carrier and non local effects, unless ad hoc phenomenological corrections are implemented.In the following, we will briefly revise the widely accepted way to deal with II in conjunction with Eqs.(1)(2)(3).

LOCAL MODEL: DESCRIPTION AND FAILURE
The most common method to deal with II henomena in Drift Diffusion simulations, is represented by the Local Model (LM) [12-15].Within such model the probability to generate a II event in (x, x + dx) for a carrier moving in the +x direction depends only on the local electric field.Under this hypothesis, II can be fully characterized by the mean free path between ionizing collisions, for both electrons, < n >, and holes, < Ip >.
The II generation rate is linked to the < li> parameters of the LM for electrons and holes respectively.Equations (4) are usually introduced referring to a constant electric field, but can be extended to a generic field shape by assuming the II coefficients to be function of the local electric field.
Within the LM, the II process in a given semicon- ductor under a specified electric field is fully charac- terized by a pair of real numbers, namely the ionization coefficients.
Equations (4-5) can be used to link the lI coeffi- cients (i.e. the inverse of the carriers mean free path between ionizing collision) to macroscopically observable quantities, such as the multiplication fac- 071 --=--Canali et al. [18]   106.
--o--Ritter et al. [17] , , i l f --e--Bulmann et al: [8].% 10 0 5 10 F " xl0 a rn V) FIGURE Dispersion of the LM ionization coefficients measured by different authors for electrons on (a) In.53Ga.47Asand (b) GaAs tor or the breakdown length [12-14].Such equations can, in turn, be applied to extract the II coefficients as a function of the electric field from experimental data [12, 13].As a matter of example, we report in Fig. la a comparison of the InGaAs electron ionization coef- ficients obtained by application of the LM from dif- ferent authors [16][17][18].The same quantity is plotted in Fig. lb for GaAs [6,8].The observation of Fig.
reveals that the ionization coefficients are dependent on the particular device under examination.Particularly, their inverse cannot be taken as a reliable estimate of the mean free path between ionizing collisions.

MONTE CARLO ANALYSIS: THE DEAD SPACE CONCEPT
The dispersion in the experimental determinations of the II coefficients shown in Fig. arises from the inadequacy of the LM.In order to clarify this point we refer to a Monte Carlo simulated experiment, by considering an electron current entering a region where a strong electric field (F=5 Position (nm) FIGURE 2 Number of ionizing collisions for unit time predicted with Monte Carlo simulation both for primary and secondary carri- ers.The quasi-exponential decay of the number of ionizing colli- sions due to the primary allows to define dead spaces and DM ionization coefficients in the simulation.In Fig. 2, the number of ionizing collisions (per unit time and volume) is plotted as a function of position.Two types of processes are con- sidered: the first ('primary') collisions of the carriers entering the high field region, and all other ('secondary') collisions, either of the primary carriers or ,of those that are generated within the field region.
Figure 2 clearly illustrates a fundamental physical feature of II process, that is the "dead-space" effect: a carrier must travel a certain distance before reaching the threshold energy for II.The "cold" injected pri- mary electrons are not immediately available for II.
Thus, a dead-space zone (dn) where no ionizing colli- sions occur is present near the contact.Equivalently, the dead-space effect can be described by the associ- ated energy (Eth,n gained from the carriers over the d n length.For a constant electric field, threshold energy and dead-space are related by the equations Eh,p q" IFl-dp; Eh,, q.IFI" d, (6) for holes and electrons respectively.It is often more convenient to speak of threshold energy because such quantity is generally less sensitive to the field value than d i, and often, as a first approximation, it is con- sidered to be constant for a given semiconductor.Eth, represents the energy that a carrier must receive from the fi.eld to appreciably initiate II.
Secondary carriers come into existence only after the primary ones collide, that is not before x d n.As it can be seen on the figure, they also cannot ionize before another dead-space, that is around x 2d n.
The reason of such behaviour can be understood in the light of the microscopy of the ionizing event.
When the primary carrier impinges on the valence electron, it release its energy almost completely to generate a hole-electron pair roughly in rest condi- tion.Accordingly, all secondary carriers can be assumed to start their motion close to zero kinetic energy.This behaviour, totally neglected in the LM, must be taken into account to adequately reproduce the II process.The mean free path between ionizing collision is only a first order description of the proc- ess.The dead space concept gives further information and, if considered in a suitable form, can add in accu- racy to an adopted model, as we will see shortly.

IMPLEMENTATION
We will refer here to a region with a step-constant, positive electric field, E We further assume, for the sake of simplicity, that carriers are approximately in the same rest condition after suffering an ionizing col- lision, or after being generated by an ionizing collision.
Under these hypothesis a very general procedure for II modeling can start by considering the spatial probability density of the ionization event, Pn(x) (Pp(x)), for a rest electron (hole) entering the field region at x 0 10].Clearly, Pn (x) 0, x > 0 and Pp(x) 0, x < 0. Field edge inclusion is easily accom- plished with the further assumption of zero II proba- bility soon after the field falls to 0. Under the previous assumption Eqs. ( 5) must be substituted with the more As we have previously shown, most of the limitations implicit in the LM arise from the fact that the dead space effect is neglected.The simplest assumption that allows us to overcome such difficulty is to consider the Pi functions to be delayed exponential func- tions of the type De -'(x-dt'), x >_ dp Pp(x) O, x < dp that is, to assume that the probability of a hole (electron) initiated II event occurring between x and x + dx (x dx) is 0 until the carrier travels a dead space length, and 5D dx (c D dx) afterwards [10].The reason for such an assumption can be supported by the obser- vation that the shape of the curve referred to as "pri- mary" in Fig. 2 is proportional to Pn(-x).Indeed, it is very close to the delayed form of Eqs. ( 9), and D and d n can be extracted from a fitting procedure.Similar arguments hold for holes.It is possible to show that the previous assump- tions, which define the Delay Model (DM) of II, lead to the following expressions for the generation rate [9] (see also [19] where carriers are used instead of currents) dJp dJn dx qG atDJn + DJp dx (10)   x-t-dn / dx' x J Jp(xdp) --jHp dx'. (12)  x-dp Equations (10-12) are valid if II is the only relevant generation-recombination mechanism.The appropri- ate boundary conditions are o, < jHn (x) --O x > W-dn Jp(O) Jpo (13) J(W) J.w The DM is actually equivalent to the model devel- oped in [10,20] to study the mean gain of avalanche photodiodes.Equations (10-12) can be easily inserted in Drift Diffusion simulators, whose solution can still be accomplished via an iterative scheme, as proposed by Gummel [21 ].Within the DM, the II characteristic of a given semi- conductor under a specified electric field are expressed by two pairs of parameters: the ionization coefficients and the threshold energies for electrons and holes, respectively.Those parameters are linked to the mean free path between ionizing through Eqs.
(7) which takes the form 1 < lp >---D + dp; < In >--D + dn. (14)  The d terms are usually negligible in the lower field range, but their contribution increases with field strength.Similar results were reported for the first time in [19] by Y. Okuto and C. R. Crowell, but up to now they have not been received considerable atten- tion.
It should be noticed that Eqs.(10) (13) have been introduced referring to constant fields, but they can be easily extended by assuming both threshold energies and II coefficients to be dependent on the local elec- tric field.

PARAMETERS EXTRACTION FOR DELAY MODEL
A practical way to obtain the DM parameters is through general techniques of error minimization.Let us suppose that N different experimental values are available, such as the multiplication factors (Mi) for electrons and holes, for a certain number of devices under a specified bias condition.First, the field shape, Fi(x corresponding to each value has to be deter- mined.Then, a sample set of the unknown functions (D(F), liD(F), Eth,n(F), Eth,p(F)) is chosen, in corre- spondence to an arbitrary number (P) of field values, F h (h (l, P)), and a gueSs is made about the range of variation of the 4.P unknown samples CD(Fh), [$D(Fh), Eth,n(Fh), Eth,p(Fh).
The extraction procedure simply consists of two further steps in the first we randomly choose the value of the 4.P unknown samples in their relative range of variation.For each choice we simulate all the N structures using a suitable function to interpolate the samples, and we calculate an error function which could be just the simple mean root square of the dif- ferences between experimental and calculated log(Mi) values.By tracking the sample set which gave the minimum value, we reduce the searching range and update its bound to maintain the central value over the best samples set.When the error function decreases under a prescribed quantity, the last step starts, con- sisting of a simple gradient search for the minimum value of the error function until the zero value is reached.
A series of improvements can be thought of in order to enhance the method convergence and accu- racy, depending on the particular device that provided the experimental data.The fundamental point is that, thanks to this method, we can successfully extract a set of highly reproducible parameters that can be used to describe II process in a given semiconductor.If we restrict our attention to structures with low overall multiplication factor (M 1 << 1), many simplifica- tions can be done on the full DM, that allow a faster parameter extraction.
We have applied the described algorithm to experi- mental data measured on npn InP/In.53Ga.47AsHBTs [17,18].In such transistors, the electron current injected from the base grows in the collector region because of II processes, the strength of multiplication depending on collector bias.The electron multiplication factor, M n, is defined as the ratio between the col- lector current and the electron current injected from the base into the collector With a suited procedure [18,22] the M n dependence on collector polarization for a given device can be extrated with great sensitivity.Both structures con- sidered in 17, 18] are operating in the range of col- lector biases that satisfy the relation M n-1 < < 1.For sake of simplicity, we have made the widely accepted assumption to consider Eth,n independent on the local electric field (at least for the range of fields considered).The procedure to extract an estimate of the OD(F) and dn(F) functions proceeds by applying the algorithm described above.We calculated the field shapes in the collector region for the device reported in 18] by means of a Drift Diffusion simula- tor, using doping data obtained with C-V measure- ments [22].For the structure reported in [17] we referred to the nominal dbping value, reported in the original work.
In Fig. 3we plot the universal electron impact ioni- zation coefficient for In.53Ga.47As extracted from our algorithm at the threshold energy of 1.25 eV in the considered field range.It is worth reminding that, if the LM was to be used (as actually done in [17,18], see Fig. 1), the ionization coefficients for the two F " xlO "a m V") In.53Ga.47Ascompared with previous LM determinations.The cal- culated values, used in DM Drift Diffusion simulator, allow a per- fect calculation of the multiplication factors in both HBTs structures would differ from one another.The calcu- lated values on Fig. 3, used in DM-based Drift Diffu- sion simulator, allow a perfect calculation of the multiplication factors in both HBTs.
In conclusion, we have shown that it is possible to extend Drift Diffusion algorithms to include non-local effects of impact ionization, thus recovering the accu- rate physical description provided by Monte Carlo simulations.
The work was partially supported by Piano Nazion- ale "Materiali Innovativi e Avanzati" of the Italian Ministry for University and Research Here o c and [L represent the ionization coefficients --Neviani et al,[6] the expected value calcula- tion.Different choices for the Pi(x) functions would lead to different degrees of accuracy of the model and for the associated computational effort.If the Pi(x) functions are taken as simple exponentials of the form Le +aL'x, x<_O Pn

FIGURE 3
FIGURE 3 Universal DM ionization coefficient for electrons on