Calibration of a One Dimensional Hydrodynamic Simulator with Monte Carlo Data

In this paper we use the code Exemplar for matching a hydrodynamic 1D, timedependent simulator and the transport coefficients obtained by the Monte Carlo simulator Damocles. This code is based on the Least Square method and it does not require any a priori knowledge about the simulator (analytical form of the equations etc.). The stationary electron flow in a one dimensional n+-n-n+ submicron silicon diode is simulated.


INTRODUCTION
Electronic transport in semiconductors can be described by hydrodynamic models (hereafter HM), obtained by taking the moments of the Boltzmann transport equation (hereafter BTE): the resulting mathematical model consists of an infinite hierarchy of Partial Differential Equations expressing balance laws for the particle number n, velocity , total energy E, deviatoric stress tensor (0"), energy flux (or heat flux h) and so on, coupled with the Poisson equation.
Recently Anile and Muscato [1] presented an Extended Hydrodynamic model where the closure of the moment hierarchy is obtained by exploiting the entropy principle: this system, which is hyperbolic, consists of 13 scalar equations in the 13 unknowns (n, , E, (ij), ) and completely determines the description of the stress and of the heat flux, at variance with the other models.Such a model has been tested successfully with Monte Carlo (hereafter MC) simulations in Silicon [2].The Left H and Side of the balance equations, called production terms, represent the average rate of change of carrier total energy Qw, momentum Qp, energy-flux Qs and stress Q</j> due to the scattering of carriers with the lattice.Usually they are approximated with ad hoe empirical formula [3] or as relaxation terms [4]: this last approxima- tion leads to a serious inconsistency with one of the fundamental principles of Linear Irreversible Thermodynamics, the Onsager Reciprocity Princi- *Corresponding author.
pie.In order to tackle this problem, we expanded the electron distribution function in Hermite polynomials around the state of global thermal equilibrium (limiting ourselves to the first two polynomials), where the electrons lie close to the bottom of the conduction band.By using a well known procedure due to Grad [5] we obtained the following equations: Qw

+ + + (3)
where To is the room temperature, E0 the lattice total energy, no a reference impurity concentration (1018cm-3).Since in this approximation the distribution function is quasi-isotropic [2] the deviatoric stress tensor /is negligible and by this method we cannot extract the corresponding production term (unless we consider higher order Hermite polynomials).
For the sake of simplicity we modeled the stress production as a relaxation term: Q<ij>-<ij> We should emphasize that: the coefficients appearing in Eqs.(1-4) are not fitting parameters but rather will be extracted by MC data obtained by the Damocles code [6], in the case of parabolic spherical band approxima- tion, in order to be consistent with the Anile and Muscato model, obtained under these restric- tions; these coefficients are not functions of the positions in the device as some authors claim.
With such coefficients we obtained a closed set of hydro equations which has been solved by an adequate numerical scheme.Since the Monte Carlo method gives a stochastic solution of the BTE, the results are noisy: how does the hydro solution change for small variations of the para- meters Tw, T, a0, al, a2, b0, bl, b2 9.In order to answer this question the simulator Exemplar, developed by one of us [7], has been used.The plan of the paper is the following: in section 2 we discuss the optimization problem; in section 3 we simulate the usual n+-n-n + diode with the Damocles code and with our Extended hydrodynamic model.We discuss the range of validity of the optimization procedure for various biases and conclusion are drawn.

THE OPTIMIZATION PROBLEM
The optimization problem is well known: defining the residual function N R() -wi(f(xi, ) yi) 2   (5) i=1 where Y is the vector of the parameters, f(xi, ) is the analytical description of the event computed in xi (modeled by the vector parameter ), Yi the real event value obtained in xi and wi an appropriate weight factor.In a mathematical language the problem of the least squares method could be expressed in the form: min{R (7) where n is the dimension of the vector 7 and 7 the solution of the system, generally, non-linear VR() -0 (7) and R" the n dimensional space of the real R.If we know the analytical equation off(x, 7), we could evaluate and resolve the system vR() 0, using one of many methods proposed (e.g.Marquardt).However, if we do not know the analytical form off(x, E), it is not possible to find the solution of the problem (6) finding the solution of the system (7).This happens when we want to determine a set of parameters to model a general device with any kind of simulator: we run the simulator obtaining the macroscopic quantities but not the derivatives of the Residual function.In this case we can resort numerical methods as the simplex method and Powell's method, which minimize the Residual function, without having to know its derivatives.
The Exemplar code runs the simulator to calculate the function f(xi, ) (i= 1,...,N) then compute and optimize the Residual function with the Simplex or Powell method (according the user choice) and the procedure restarts until a tolerance of the Residual function is obtained.
The proposed algorithm can be used with any simulators which can read their input from files and store their output in column-formatted ASCII files.It was implemented in a software program with a friendly user interface based on X-Window, which resolves the problem of the optimization well, but, since during the optimization it runs the simulator, it needs a large computational time.The n + -n-n + diode consists of two n + regions 0.1 gm-long doped to a density ofN 1018 cm-3, while the central n region is 0.3 gm wide, with a doping density of N 1016 cm-3.We consider Ohmic boundary conditions, To 300 K, and V of applied bias; simulations refer to the stationary regime.
In Figures and 2 we plot the MC data for Qv and Qs and the fitting with eqs.( 2) and (3) (with ooo): we see that our functional form fits well the data.We observe that the MC data are noisy especially near the contacts (this run needed month of CPU in a IBM Risc 6000 590) and consequently the fitting coefficients are not unique.By using these coefficients the Extended Hydro- dynamic model is solved by using a simulator based on the splitting method between relaxation and convection (Tadmor scheme) [8].Then we run the Exemplar code in order to find the best coefficients which fit well the MC data.In Figures 3, 4 and 5 we compare the MC data (with ***), the hydro data without optimization (with ooo), the hydro data optimized with Exemplar (with xxx).The optimized coefficients (shown in Tab.I) differ from the previous ones at maximum by 40%.We notice that the well known 'spurious' peak in the velocity curve Figure 3 is lowered and the energy and heat flux curves Figures 4 and 5 are closer to the MC data.However this 'spurious' peak cannot be eliminated by any small change in the parameters and therefore its persistence calls for a more radical examination of the basic assumptions of the model.
A crucial question to be addressed is what is the range of validity of the optimization procedure?In order to answer this question we simulate the device with the coefficients of Table I (obtained with Volt) for various biases.The result is that in the range 0.8 / 1.2 Volt the hydro simulations fit well the respective MC data (see Figs. 6, 7 and 8 in case of 1.2V): for higher biases the hydro simulations diverge considerably with respect to the MC data.Finally we say that the modeling of the production terms and the closure of the hydro model are crucial points.They should be obtained by using the physics (e.g., an entropy principle, expansions of the distribution function).The transport coefficients can be obtained either by experiments or by MC simulations, and their validity is restricted to a neighborhood.
We are trying to improve our model and to simulate 2D devices: work along this line is in progress and will presented elsewhere.

FIGURE 2
FIGUREAverage rate of change of carrier momentum Qp vs. distance: Monte Carlo data (with ***) and the fitting formula Eq. (2) (with ooo).

5 FIGURE 3 FIGURE 5 FIGURE 4
FIGURE3 Average drift velocity vs. distance computed by using the hydrodynamic simulator without optimization (ooo), the hydro data optimized with the Exemplar code (xxx) and Monte Carlo data (***), with Volt bias.

TABLE
Transport coefficients (in psec-), channel 0.3 tm, Volt bias Average electron energy vs. distance computed by using the hydrodynamic simulator with transport coefficients of Table (obtained with Volt) (ooo) and Monte Carlo data (***), with 1.2 Volt bias.Heat flux vs. distance computed by using the hydrodynamic simulator with transport coefficients of Table (obtained with Volt) (ooo) and Monte Carlo data (***), with 1.2 Volt bias.