Applicability of the High Field Model: An Analytical Study via Asymptotic Parameters Defining Domain Decomposition

In this paper, we present a mesoscopic-macroscopic model of self-consistent charge transport. It is based upon an asymptotic expansion of solutions of the Boltzmann Transport Equation (BTE). We identify three dimensionless parameters from the BTE. These parameters are, respectively, the quotient of reference scales for drift and thermal velocities, the scaled mean free path, and the scaled Debye length. Such parameters induce domain dependent macroscopic approximations. Particular focus is placed upon the so-called high field model, defined by the regime where drift velocity dominates thermal velocity. This model incorporates kinetic transition layers, linking mesoscopic to macroscopic states. Reference scalings are defined by the background doping levels and distinct, experimentally measured mobility expressions, as well as locally determined ranges for the electric fields. The mobilities reflect a coarse substitute for reference scales of scattering mechanisms. See [9] for elaboration.


INTRODUCTION
This presentation is motivated by the search for more reliable macroscopic models of high field transport in submicron structures which are computationally efficient.This is a research topic of major interest to the micro-electronics industry.
Transport in submicron structures differs from transport in bulk material in many ways.It includes far from equilibrium situations created by large electric fields in small structures; ballistic electrons, when dimensions of the drift-forcing term are of the same order as the mean free path; and new effects, introduced by the close proximity to the boundaries of regions where the density gradients are (relatively) rapidly changing.
Most of the previous work on transport in submicron semiconducting structures used several corrections to drift-diffusion equations in which the current is a sum of drift and diffusion terms, where both have field-dependent mobility and diffusion coefficients, respectively.References on models and computational methods for N +-N N + structures of Si or GaAs can be found in the literature of the last decades [24, 23, 16, 21, 25, 20].However, the computations have not yielded results which have been genuinely satisfactory, and the main reason has been the lack of under- stading of the accuracy of the modeling in many heterogeneous structures.An alternative to the drift-diffusion equations has been Monte Carlo simulation, a very costly computational procedure.
An approach intermediate between drift-diffu- sion and Monte Carlo simulation was generated by the hydrodynamic like theory, using velocity mo- ments of the Boltzmann equation (see [6, 2, 3, 15]) and realizing closure by imposing constitutive relationships between the state macroscopic vari- ables.Extensive numerical modeling of N + -N-N + Si and GaAs structures under several geometric constraints can be found in [12, 17, 18, 19].For the derivation of these models, the distribution function is often assumed to be a Maxwellian distribution.However, this assumption has been shown to be questionable in lightly doped submicron structures, in part because of the presence of ballistic electrons.
However, it has been mathematically proved that hydrodynamic models relax to energy trans- port models.The latter have been derived as macroscopic limits associated with a particular choice of space-time scale that makes some collision mechanisms dominant (see [10] for a rigorous proof in the case of the relaxation of the full hydrodynamic model).These models, how- ever, cannot describe situations where the drift- forcing term and the resulting scaled mean free path are of the same order, but they are asymptotically correct in regions where the predefined background is very slowly varying (i.e., almost constant doping).A new approach at- tempting to fill in the middle ground among drift- diffusion, hydrodynamic energy transport, and Monte Carlo simulation (see [4, 5, 26]) was the use of direct, numerical solution of the Boltzmann equation within the relaxation time approximation.The use of the Boltzmann equation for non-electron transport is justified in the above references, as well as in [7,13].The relaxation time approximation within the Boltzmann-Poisson sys- tem approach, positing that all scattering pro- cesses can be characterized by a few scattering rates setting up the scale of relaxation time constants, may not be strictly valid for GaAs [11].
Our goal, however, is to present a full solvable model for an N+-N-N + structure, as an alternative to all previous calculations, which will produce comparable results with greater efficiency than attained by [4,5,20].In particular, we incorporate into the modeling the macroscopic limit associated with the choice of space, time, potential drop, and local electric field scales, combined with the relaxation time scale, which makes the drift-forcing term and collisions both dominant, and of the same scale set up by the local scaled mean free path.
As a consequence of this approach, we incorpo- rate a level of modeling that takes into account mesoscopic-macroscopic multiscales, and corre- sponding limit derivations, according to different scales becoming dominant in different regions of the N + -N-N + device, depending on the spatial inhomogeneity of the lightly doped region.

PRELIMINARIES OF THE HIGH FIELD MODEL
We identify three dimensionless parameters, viz., The ratio r/of drift and free velocity (the latter usually taking on the thermal velocity); . .The scaled mean free path e and The scaled Debye length -,/from the electrostatic potential equation of a self-consistent model.
The dimensionless parameters induce domain dependent macroscopic approximations, which are valid within regimes that do not change the scale of the Debye length.This is of paramount importance in the macroscopic derivation asso- ciated with the BTE under a self-consistent field: the scaling of the forcing term depends on the scaled external and internal fields, and must not override the limiting regime of the macroscopic derivation.Therefore, this model incorporates kinetic transition layers linking mesoscopic and macroscopic states.
Reference scalings are defined by the back- ground doping levels and distinct, experimentally measured mobilities which serve to set up the scale of the relaxation time approximation of the scattering mechanisms.In addition, scales take into account locally determined ranges for the electric fields.Numerical experiments of this proposed model are also presented by the authors in [8].Although this is a first approach to domain decomposition models, based on a mesoscopic- macroscopic linking, our results seem to give an improvement (see also [13]) of computational efficiency, relative to the kinetic computations of Baranger and Wilkins [5].Moreover, the model itself seems to be an improvement over the augmented drift-diffusion model used by Kan, Ravaioli and Kerkhoven [20].This model, based on a model introduced empirically by Thornber [25], employs necessary corrections to the usual macroscopic approach to high field models.

BOLTZMANN-POISSON SYSTEM SCALES
Finally, we concentrate on the derivation of the macroscopic models for semiconductors, based on asymptotic expansions of multiscale dimensionless BTE-Poisson systems.For the sake of simplicity we shall consider only electrons.The semiclassical Boltzmann-Poisson system, within a parabolic band-relaxation time approximation, may be written: where F(x, v, t) is the density function for the electron at position x, velocity v and time t.The constants e, m,and e represent the electric charge constant, its mass, the relaxation time reference scale, and permittivity of the material, respectively.The functions b and E of (x, t) represent the electrostatic potential and its field.
The bracket (.) denotes the usual average of the distribution function with respect to the velocity variable.Finally M= M(v, x) denotes the Max- wellian centered at zero velocity with 0 O(x) representing the background temperature of the lattice, measured in units of specific energy.Derivations of mesoscopic limits associated with system (1) depend on a particular choice of space- time scales, collision mechanism scales, and drift velocity scales induced by self-referencing field scales.These regimes, as we shall see, can make some terms of the BTE dominant.The idea that follows then, is to expand the solution of the BTE about the solution of the grouml equation corre- sponding to the dominant terms.Of course, this will depend on the form of scattering mechanisms that set up the scales of the collision terms.Lately, a very nice survey of a hierarchy of macroscopic models for semiconductors under limit regimes that correspond to dominant collision mechanisms has been presented in [1].These macroscopic models correspond to classical drift-diffusion and so-called energy transport models.They arise from expansion about kernels of the collision operators, which are field independent, so they do not incorporate the scales of the external, or even self-induced electric field.Diffusions are field independent, and consequently they cannot ap- proximate the ballistic electron distribution func- tions that have been observed by Barenger and Wilkins [4,5].After the work initiated by [14] and [22] on strong forcing scaling for external fields, where the dominant term in the scaled Eq. ( 1) is given by the balance of the force-drift (elm) E(x, t)X7vF and the collision terms, the following 3- scale dimensionless formulation has been formu- lated in [9] in g Euclidean dimensions: rift + v. VxF +rlvxO VvF--{(F) M-F} zxO (e) Ue( ) (2) 7"01/2 is the scaled mean free path for a where e T length scale L and relaxation time r, and 0 is the reference scale of the background energy.The dimensionless constant r-U/O 1/2, where U, in 7"e, takes into units of velocity, is given by account the reference scale of the local electric eL field.Finally, the dimensionless constant -y-,.[4] represents the scaled Debye length, where scales the density of the fixed background Nd(x).Now, a strong force regime corresponds to values of r/ O(1), -y O( 1) and e < L. In this case, system (2) is obtained, where the distribution function F is scaled with the drift velocity U (instead of the low field scalings that use the thermal velocity scale for F).The time scale is then fixed at LU -.

ASYMPTOTIC LIMITS AND THE HIGH FIELD MODEL
Clearly, any asymptotic expansion is justified under the assumption that the scaling does not break down.Thus, any expansion of the BTE-Poisson system (2) requires that the scale of the internal field defined by the Poisson equation (i.e., depending on the value of 3') does not override the scaling of the BTE. is the reference scale of the predefined Na(x) and will change scale whenever Nd(x) does.Accordingly, so will the local field.Therefore, as an immediate consequence of the fact that -y q,(#) can be easily checked for GaAs (as in [5, 20]), the dimensionless quantities vary in the N+regions from the N regions.Thus, in N + re- gions, r/ e o(1).Because of scaling times and the distribution function F, accordingly, (2) yields the Drift-Diffusion-Energy Transport regime as the relative electric field remains moderate and the collision mechanisms become dominant.
However, in the N region it is computed by the authors in [8] that e (1/7) and r/ O(1).Thus, while 3' remains of order O(1), it has been shown that densities and the electric field satisfy the equations, pt-t-V J-O, J---T[Vx(lg(O q-2E @ E))] #oE e (#,Uw), +'r#EpV.(#e) "r#p-27.w O, curl w #E (R) 7p, E --VO, (3) in their dimensionalized formulation.Here # is the mobility.The new variable w can be interpreted as the curl of a magnetic field, associated with the high electric current.This model is a correction of the one proposed by Thornber  [25], where the diffusion and transport coefficients were given empirically.
We point out that in this relaxation time approximation, the mobility coefficient is given by # , so clearly the saturation of the mobility depends on how the scattering rates scale for r changing with respect to the electric field.The standard assumption based on curve fitting is to take #-C/(1 + C21EI2) 1/2 (see [16]).
We stress that this regime is valid only for the strong force scaling.The transition regime from weak to strong forcing scaling, and conversely, ought to be done by carefully considering the corresponding kinetic transition layers.The im- plementation and solution of this transition problem require solution of boundary value problems at the kinetic and macroscopic level as well.

DOMAIN DECOMPOSITION AND TRANSITION LAYERS
Appropriate boundary data at the kinetic level and the corresponding limiting data for the mesoscopic approximation are presented in [9].For kinetic level distribution functions there are prescribed exact neutral space charge conditions at contact boundaries and standard insulating boundary conditions at insulating walls.Then, these condi- tions, together with the asymptotic approxima- tions, yield Dirichlet boundary conditions at the macroscopic level, with a correction term of the order of the asymptotic parameter (i.e., the order of the scaled mean free path e).The contact condition, at the kinetic level, on a space boundary section I" reads [9]: where n denotes the outer normal with respect to the boundary section F, and P denotes the lowest order term in the expansion for either weak force scaling (i.e., P is a Maxwellian) or strong force scaling (i.e., P peaks asymmetrically and is electric field dependent).
The corresponding macroscopic condition for the inner boundary can be computed explicitly [9] and yields Here, p(x) Na(x) + o(e)K(x), We finally remark that the Ohmic contacts on N+regions correspond to a local scaled mean free path e of much smaller order than in the N regions.Hence, at boundary ohmic contacts on N + regions, condition (5) is almost a neutral space charge condition at the macroscopic level, as historically has been, and still is typically pre- scribed.However, at the transition layer from the N + to the N regime, the transfer of the data entails a kinetic computation using condition (4); then, one prescribes an exact background density, and thus delivers, at the next level, for the strong force model, a condition for the density that is an O(e) derivation from the background density.This fact can actually be observed in the computations of Barenger and Wilkins [4,5].Numerical implementations of strong-weak forcing decomposition are presented by the authors [8] in this same issue.However, we use there a matching with hydrodynamic computations of well accepted performance.In Figure 3 of [8], we present a graph of r/, which clearly correlates with the suggested domain decomposi- tion.Theoretical justifications and alternative numerical methods of these domain decomposition ideas are underway and will be presented in future work.