Approximation of the BTE by a relaxation-time operator: simulations for a 50nm-channel Si diode

In this work we present comparisons between DSMC simulations of the full BTE and deterministic simulations of a relaxation-time approximation for a nowadays size Si diode. We assume a field dependent relaxation time fitted to give the same drift speed (mean velocity) as DSMC simulations for bulk Si. We compute the density, mean velocity, force field, potential drop, energy and I-V curves of both models and plot the pdf of the deterministic relaxation-time model. We also compare the results to augmented drift-diffusion models proposed in the literature to approximate the relaxation time system in the quasi-ballistic regime. The quasi-ballistic and ballistic regimes are distinguished by using local dimensionless parameters.


INTRODUCTION
In the present report we consider an electron-hole transport model for a 250nm n+-n-n + Si diode with a 50 nm channel.Due to the small value of the doping of holes, we assume a constant doping of holes over the whole device and we compute only the evolution of the electrons.This assumption simplifies the computations and they are qualitatively equivalent to the more realistic double carrier model.We refer to [1] for more details about this simplified and idealized device.One of the goals of our studies is to establish rigorous comparisons between probabilistic computational methods (DSMC) versus deterministic ones for the BTE in the modeling of nanodevices.Our direct simulation of the BTE involves a selfconsistent field dependent relaxation approximation.In addition we compare with our Drift-Diffusion-Poisson (DDP) solver and a recently implemented domain decomposition technique to compute non-equilibrium statistics in the channel region.We use a standard DSMC code based on classical texts [7].For studying this evolution, we consider the full BTE Of # 0-Vfx --m E(t, X)fv Q(f) in the parabolic band approximation.Here, f= f(t, x, v, w) is the density function for an electron at position x E [0, L] and velocity (v, w)E R 3 (v is the first coordinate in velocity) at time > 0, where L is the device length.The constants e and m represent the unit charge and effective electron mass, respectively.The electric field E= E(t,x) is self-consistently produced by the electrons moving in a fixed background with density of holes (acceptors) na(X) and electrons (donors) rid(x).
Thus, E is determined by the Poisson equation where now f depends only on (t,x, v) and the density p and the current j are computed by taking integrals only over v.The distribution Moo is the classical absolute Maxwellian centered at zero velocity with variance given by Oo=(ke/m)To, with the Boltzmann constant ke and the lattice temperature To in Kelvin.Our goal is to make a comparison between the results of the full 3-D velocity BTE solved by a DSMC method [7] with the simulation of the 1-D relaxation-time equation (1.2) solved by deterministic WENO methods [4,6].
The relaxation timedepends on the absolute value of the force field in such a way that the drift speed #(E)IE approximates the one computed by the DSMC solver.Thus, #(E)[E] =(e/m)-(E)lE is linear for small values of IE] with slope #0 and saturates at Vd as ]E becomes large (see Fig. 1).
We use the classical constitutive relationship eoxx e(p(f) nd(x) 4-na(X)), E(t, x) -'x, e 2#0 1) where eo is the permitivity of the material and p(t, x) f f t, x, v, w)dvaw, J j(t,x) /3 vf(t,x, w dvdw x e [O,L], t_>0 are, respectively, the charge and the current (in the field direction) of the electrons.We will consider acoustic and non-polar optical phonon scatterings which are the most relevant scattering mechanisms in Si [7].
The classical kinetic relaxation model for charged transport is derived from a low density approximation of the semi classical Boltzmann-Poisson system.Though this system is posed in the three dimensional velocity space, a one dimen- sional model recovers the important features of the charge transport that are given in the direction parallel to the force field.Therefore this equation reads [5] Of The parameters #0 (low-field mobility) and va (saturation speed) are adjusted calibrating the drift speed in the bulk Si with respect to the Monte Carlo code.Accordingly, we take #0=0.19 and va=0.15.In addition, DSMC SIMULATION 351 compared to (1.3).In this way, we make sure both models refer to the same physical device.Our simulations of the field dependent relaxation BTE, strongly indicate that most of the channel device operates in a quasi-ballistic (or Drift- Collision balance) regime where the drift transport term balances the collisions.For these transport scales, the field dependent relaxation-Poisson system (1.2), (1.1), and (1.3) can be well approxi- mated by drift-diffusion like models inside the channel producing faster numerical solvers for non-equilibrium regimes.As it was worked out in [2], one can use a hybrid domain decomposition technique based on the scales involved.This tech- nique consists of the computation of classical drift-diffusion-Poisson (DDP) models in the n + regions which correspond to dominant collision processes, and augmented-drift-diffusion-Poisson (ADDP) models in the n-channel region.The procedure to choose the hydrodynamic set of equations, used for approximating the 1-D relaxa- tion time system, is based on the values of two dimensionless parameters c and /. a--B0' /3--B0 (1.4) where Uo and Bo are the local drift and ballistic speed, respectively.Assuming that a /3 0(1) (quasi-ballistic regime), one can derive [3] the (ADDP) system.The ADDP system consists of the continuity equation where the current is approximated by ( e )P(-#PE + w) coupled with (1.1) where w (#pE)Ix=x and x is some point in the computational region.The total current j(t,x) is approximated by the J(t,x) solution of (1.5).We refer to [3] and the references therein for a deeper discussion of this system and the role of w.The density p(t, x) and current J(t, x) that solve the hybrid DDP and ADDP system are asymptotic approximations (Chapman-Enskog) of the first and second moments of f(t,x, v) solution of the kinetic system (1.1)-(1.3) in the n+ and n- channel regions respectively.
However, if c=0(1) and /30 (ballistic re- gime), then the derivation of a hydrodynamic set of equations is done by a closure procedure [2].Let us remark that in this particular device the ballistic scaling does not appear due mainly to the small value of the saturation speed.

NUMERICAL SIMULATIONS
We consider a one dimensional Si n+-n-n + structure of length 250nm with a n-channel of 50nm.Here, the donor doping profile given by rid(X) is a sharp step function with density values 5 10:Z4/m 3 in 0 <_ x _< 100 nm and in 150 < x <_ 250; and 1021/m 3 in 100 <x < 150.Also, we consider a doping of holes constant over the whole device of value na(X)-5 1022/m 3.
All the simulations for the deterministic models are by the fifth order WENO finite difference schemes [4,6].These are explicit schemes which are extremely stable and robust, especially suitable for strong convective regimes.

CONCLUSIONS
When comparing IV-curves of a 50 nm channel Si diode (Fig. 2), the deterministic simulation of a calibrated field dependent-relaxation compares very well to DSMC simulations up to bias of 0.1 V.For larger bias the relaxation BTE simula- tion fails to capture the stronger scattering effects  Normalized pdf of the relaxation 1-D system (solid line) for the Si n +-n-n + device at Vbias--0.6V at 132 nm compared to the dominant pdf P (dotted line) of the quasi-ballistic regime.x FIGURE 4 Numerical comparison of density, mean velocity, force field and potential and IV-curves for the deterministic 1-D BTE, DDP and ADDP simulations for the Sin +-n-n + device at Vbias 0.6V.Solid line: the relaxation 1-D BTE; dotted line: the ADDP system, squares symbols: the DDP system.
modeled by the whole collision operator as well as the 3-D effects.This effect is well illustrated by Figure 2, where the density, mean velocity, force field, potential and energy for both solvers have been plotted for 0.6V bias.It is clear the relaxation-BTE model grossly underestimates the mean velocity and the energy.
The testing of the hybridization of the domain decomposition method of [2] produces a very good agreement for the IV-curves (Fig. 2) with respect to the relaxation-time and improves Drift-diffu- sion which clearly overshoots the DSMC and relaxation simulations for all positive bias.We use the ADDP system to approximate the relaxation- time system because the dimensionless parameters are almost equal over the channel (Fig. 3).
The hybrid-domain decomposition implementation (Fig. 4) improves the density, force field and potential with respect to the DDP results; how- ever, the mean velocity is only slightly improved.
In fact plots of the pdf in Figure 3 show that the quasi-ballistic regime is well achieved at the channel location x0--132 nm.The formula for P can be found in [2].This is a strong indication of an asymptotic non-equilibrium fluid regime approximating the kinetic relaxation-BTE.
It is important to notice that we are comparing 24 hour running time for a DSMC simulation versus a 45 minute for a WENO simulation of the 1-D BTE and a 3.5 minute computation of the hybrid domain decomposition method (DDP-ADDP).
However it is quite surprising that, in appear- ance and without asymptotic justification, the Drift-Diffusion IV-curves simulation appear closer than that from the ADDP to the DSMC one for larger bias.The reason for this is that the ADDP system was obtained as asymptotics for the relaxation-time operator.A similar approach based on the full BTE might produce a different corrected drift-diffusion-like system able to com- pare better to the full BTE.

FIGURE
FIGURE Drift speed for bulk Si.Solid line: #(E)E[ according to Formula (1.3); dotted line: the full BTE by DSMC.