A Spectral Flux Method for Solving the Boltzmann Equation

A spectral method for solving the Boltzmann equation by the scattering matrix approach is presented. The algorithm discussed can be used to simulate both bulk and device properties with arbitrary field profiles. Although the primary goal is to reduce the data storage problem of the scattering matrix approach, many of the concepts and mathematical properties developed may be useful for other traditional spectral methods as well.


I. INTRODUCTION
Nonequilibrium effects such as velocity overshoot and hot-electron injection over barriers are becoming increasingly important for submicrometer devices.' The traditional drift diffusion equation, which is based on a quasiequilibrium approximation of the Boltzmann equation, does not simulate these hot-electron properties well.During the last few years, therefore, a number of different schemes have been proposed to solve the space-dependent Boltzmann equation for high-field regions.These schemes include the Monte Carlo tec.hnique,' spectral methods,3>4 hydrodynamic approaches,s'6 etc.Each of these techniques has its own merits and limitations; for example, the Monte Carlo method is remarkably accurate, but it is computationally inefficient and inherently noisy; the hydrodynamic approach is very efficient, but it does not have accuracy comparable to Monte Carlo.Spectral methods provide direct solution to the Boltzmann equation and do not make as many untested assumptions as hydrodynamic codes do.So far, a number of promising results have been reported3 using this approach; however, this technique has not yet found wide use for general purpose device simulation.
Recently, a new technique called the scattering matrix approach (SMA) has been proposed.7'8 It shows accuracy comparable to Monte Carlo while retaining the elegant flexibility of the hydrodynamic codes.The modular nature of this approach makes it suitable for semiclassical as well as quantum-device simulation.However, the SMA in its present form has a major limitation: It requires considerable computer memory.In this approach, one precomputes a set of scattering matrices for a thin semiconductor slab under different electric-field strengths and stores the matrices as a library to be used subsequently in device simulation.The elements of these scattering matrices represent the scattering of particles from one momentum st.ate to another while transmitting across the slab.For realistic device simulation, momentum space must be resolved into a large number of rectangular bins.Since the number of scattering matrix elements increases as the square of the number of momentum bins, storing them in computer memory becomes increasingly difficult.To reduce the library size, therefore, it is desirable to have a more effic.ientscheme of momentum space representation in terms of basis functions.
The choice of the basis function to represent the momentum space is motivated by the shape of the distribution function.gWith a proper choice of basis functions, one should be able to represent the distribution function with relatively few coefficients.The basis functions available for simulation purposes include numerical basis functions, nonorthogonal basis functions, orthogonal basis functions, etc.The physical motivation of using orthogonal polynomials as basis functions is that at low fields the distribution function is approximately Maxwellian, which can be adequately represented by a few Hermite polynomials.Also, the coefficients of the orthonormal functions are inherently optimized in the least-squares sense.This is the essence of the spectral method.At higher fields, the distribution function becomes highly asymmetric and one needs an increasingly greater number of coefficients to represent the distribution function.
The purpose of this article is to investigate the issues involved in a spectral method for fluxes and to study its usefulness.We show that the required algorithm is not trivial, indeed, the basic ideas presented here may prove useful for other spectral approaches as well.The rationale behind the choice of the basis function and the completeness of the chosen set is discussed in Sec.II, along with a discussion on the orthogonal transformation.In Sec.III we show that even though a simple c.hange of basis functions is suffic.ientto simulate bulk characteristics, the transformed matrices do not obey the usual cascading rules; therefore, the matrices must be properly transformed for device simulation.In Sec.IV we discuss the computation of various quantities of interest such as carrier concentration, average velocity, etc.In Sec.V, we present some example calculations and in Sec.VI we summarize and conclude with a brief discussion.where J' and J-are particle fluxes represented by PX 1 column vectors resolved into bins in momentum space and t+, r; ) t-, and r+ are PXP submatrices with elements describing the particle exchange among various bins of momentum space.Here P is the number of modes in momentum space and the electric field is assumed to be oriented along the negative z direction and electron transport is The desired basis functions for the SMA should have assumed.By definition, Jt is the positive flux stream with two properties: (i) The functions must be defined over half of momentum space along the k, direction, since the funo tions to be expanded are flux functions; and (ii) the basis k,> 0 and J-is the negative flux stream with k, < 0. Therefunctions should have basically a Maxwellian distribution fore, the tlux function is defined over the semi-infinite space over their independent variables kZ and k; however, none of the orthogonal basis functions arising from Sturm-of momentum space in the z direction.This property of the Liouville equations simultaneously satisfy both those requirements.For example, Laguerre polynomials are de-flux functions will make necessary a few modifications of fined over the appropriate interval, however, they do not have the correct distribution over momentum space.Her-the standard spectral method discussed in the literature.3mite polynomials, while having the correct distribution function, are orthogonal over the interval 00 to -to.
To resolve this problem, we fictitiously extend the positive (negative) flux function to the negative (positive) velocity range making it an even function in velocity.This redefined flux function can be expanded by a basis set consisting of only even-order Hermite polynomials (EOHP).Bach element of this set is a product of two even Hermite polynomials, one along the longitudinal momentum k, and the other along the transverse momentum kp The set is both orthogonal and complete for all possible flux func-tions.A proof of these properties is given in Appendix A. The odd basis functions do not contribute anything to this expansion because their overlap integral with the even flux functions are zero.

Orthogonal transformation
Let us simplify the notation of Eq. ( 1) by labeling the scattered fluxes as J, the incident fluxes as Jin, and the scattering matrix as [Ml.In this new notation, Eq. ( 1) becomes (4) = EM 1 (Jin).In principle, the matrix [B] is an infinite matrix because the elements of the set of EOHP are infinite.In practice, however, one can do away with all but a few polynomials.Thus the matrix [BJ turns out to be a rectangular matrix of order NXP, where N is the number of polynomials and P is the number of rectangular grid points.Since N can be a very small number, one can now store a much smaller matrix [N'], reducing the computer memory problem discussed in Sec.I.However, since the basis function is no longer rectangular, we shall face some unique problems in its use for device simulation, and straightforward application of the algorithm present in Ref. 7 is no longer possible.We address this issue in the following section.

Iii. SIMULATION PROCEDURE
A. Bulk simulation Bulk material properties were simulated in Ref. 7 by using periodic boundary conditions.Mathematically, this is equivalent to fmding the eigenvector of the scattering matrix corresponding to the eigenvalue 1.' Therefore, so-lution of the bulk properties from reduced library matrices entiate it from the orthogonal coefficient matrix (OCM) [Jf'] involves two steps.First, the matrix is solved for [.'I.The OFM is Markovian and preserves flux according eigenvector corresponding to the eigenvalue 1, i.e., to the standard cascading rules.
[M'] (JL, = (Jk,. Second, one returns to a rectangular basis function using the transformation ub)=[BIu& (5) where J, is the eigenvector in the rectangular basis and JL is the eigenvector solved from the reduced matrix l&f'] itself.One can compute all relevant bulk properties from Jb using the Eqs.( 8)-( 12) of Ref. 7. One should note that step 2 of the procedure is not necessary, because from the vector Ji one can in principle solve for all the physical parameters.However, use of step 2 simplifies calculation and evaluation of various integrals.
It should be noted that if the-number of polynomials is small, then the column sum of [M] will not be exactly 1 as required by the cascading rules of Markov matrices.The reason is the numerical error associated with describing the emerging fluxes by a small number of orthogonal polynomials.In these cases, a proper resealing of the coefficients will be sufficient for flux convergence.
Once iteration is completed, a vector is specified in between each pair of semiconductor slabs.In order to relate the coefficients of these vectors to the fluxes in between the slabs, we need to establish the following matrix property.Let Fj and c> be the elements of the incident flux vector for OFM and OCM, respectively.These elements are related to each other by B. Device simulation For a complete solution of the device properties using the SMA, the device is first divided into many thin slabs and then a constant electric field is assigned to each of these slabs in such a way that it approximates the field profile of the entire device.The scattering matrices corresponding to these thin slabs are cascaded to simulate steady-state properties of the device.The rules for cascading matrices are given by Eq. ( 2 We shall the follow the second approach.
In this subsection, we shall state the key results only.Details of the schem_e are given in Appendix B. First, we transform [M'] to [LM] by where the matrix elements are given by where m;i is the transmission coefficient from the orthogonal modej to the orthogonal mode i and Wi and 1~1~ are the areas under the curve of the ith andjth elements of EOHP, respectively, i.e., --cc --m The new scattering matrix [$I given by Eq. ( 7) will be referred to as an orthogonal flux matrix (OFM) to differ-Fj=C>(Wj/Wl). ( 9) Corresponding to these incident fluxes, the resulting scattered fluxes Fi and q[ are given by c 10) and ( 11) Using Eqs. ( 7) and ( 9), one can show that the elements of the scattered fluxes, j$ and qi are related to each other by Fl=q;(uvwl).( 12) This important relation helps to translate the results computed using one set of basis functions to those corresponding to the other basis set.A corollary to this pr0pert.y is that the eigenvectors corresponding to eigenvalue 1 for the two matrices given by Eq. ( 7) are related to each other by F+=Ci(W/Wlj, ( where & and ci are the elements of the eigenvectors of the orthogonal flux matrix and the orthogonal coefficient matrix, respectively.Once the elements of the vector JL, i.e., q;, are obtained by repeated application of the above properties, one can use Eq. ( 5) to compute fluxes in rectangular basis function and subsequently obtain all relevant physical parameters using Eqs.( 8)-( 12) of Ref. 7. The spectral flux method can be briefly summarized as follows.We begin with a 2Px2P scattering matrix [M] computed for P rectangular momentum bins.A set of N'gP even-order Hermite polynomials is then selected, and a ~Mx 2N coefficient matrix [M'] is evaluated from Eq. (3 ).Since [1M'J is not Markovian, the normal cascading rules for the scattering matrices do not apply, so [M'] is transformed to [M] according to Eq. ( 7).The transformed Markovian scattering matrices are then cascaded, and the steady-state fluxes v_ersus position are obtained.From the steady-state fluxes Jb, we find the coefficient vector from Eq. ( 12).

IV. RESULTS
To illustrate the spectral scattering matrix approach, we present some sample calculations.First, we consider simulation results for bulk silicon.In Fig. 2, we show the velocity versus field curve computed from the spectral method.We also show Monte Carlo results based on Ref. 11.Three and five coefficient spectral methods show reasonable agreement with the Monte Carlo simulation.For the energy versus field curve shown in Fig. 3, agreement with Monte Carlo data is less satisfactory.However, at low fields the agreement is much closer.At low fields, the distribution function is almost symmetric, therefore a few coefficients are sufficient to represent the distribution function.However, at high fields one needs more coefficients to maintain the y;lme level of accuracy because the distribution function becomes more asymmetric.A comparison of the five-coefficient energy versus field curve to that of the three-coeilicient energy illustrates this point.
Next we present a simple non-self-consistent simulation of high-field electron transport for a model silicon device whose field profile is shown in Fig. 4(a).In Fig. 4(b), the velocity profile is shown as a function of position.If we compare the results with Monte Carlo simulation, we see that nonstationary transport in the device is well simulated.However, the five-coefficient spectral method is not suitable for energy simulation at very high fields.More coefficients will systematically improve the agreement.V. CONCLUSIONS We have presented a spectral flux method for solving the Boltzmann equation within the SMA framework.This method was shown to be useful both for space-dependent and space-independent simulations.Systematic improvement of the simulation results is possible by increasing the number of coefficients of the orthogonal polynomials.
For the scattering matrix approach, the significant reduction of the memory size is possible for fields below 10 kV/cm.At higher fields, due to the asymmetry of the distribution function, the saving in memory may not be significant.
Instead of using orthogonal basis functions, as we have done in this article, it is possible to use numerical nonorthogonal basis functions for device simulation.One possible choice could be a set consisting of bulk flux functions computed for different electric fields.Fluxes corresponding to an intermediate field can be obtained by a linear interpolation of the elements of the numerical bulk flux functions.The only problem of this approach is that since the numerical basis functions are not orthogonal, the expansion coefficients will not be optimized in the leastsquares sense.lo Optimizing the coefficients in the leastsquares sense will add additional complication to the process.However, once the coefficients are obtained, the device simulation algorithm developed in this article will apply regardless of the choice of the basis functions.Therefore, we have chosen to concentrate on orthogonal polynomials to clarify the basic concepts involved in using the spectral flux method.Note that Eq. (Bl 1) follows from Bq. (BlO) by using Eq. ( 8) and by using the definition of [M'] from Eq. (3).Although, we considered only the transmission submatrix for illustrative purposes, such 3 relation can be proved for all other submatrices of the scattering matrix of Eq. ( 1) 3s well.This, therefore, completes the proof of Eq. ( 7).
FIG. 1.A schemstic diagram of the incident and scattered fluxes of a thin semiconductor 4nb.[q represents the scattering matrix that relates the fluxes.
Ref. 7, a rectangular basis function was assumed and the elements of [M] were computed using Monte Carlo simulation by injecting a large number of particles within each bin of momentum space and then counting the number of transmitted and reflected particles in other bins.Let us also define an orthogonal matrix [B] with each column defining one element of the EOHP set.The required orthogonal, similarity transformation is given by where Ji represents the coefficient vector of the scattered fluxes in Hermite polynomial representation and [M'] denotes the corresponding coefficient matrix.One should note that an orthogonal transformation preserves the trace and the determinant of a matrix because the eigenvalues are not changed under such a transformation.'*This property will be very useful in our subsequent analysis.Note that the orthogonal transformation has an elegant physical interpretation.One can evaluate the elements of [&?] by injecting a flux of carriers weighted by the normalized orthogonal basis function into a thin semiconductor slab and resolving the outscattered fluxes (both transmitted and reflected) into elements of the EOHP set.
) in Ref. 7; however, these rules apply only to Markov matrices as discussed in Ref. 8. Since the orthogonal transform of the Markov matrix [M] is not Markovian (its columns do nut sum to unity), these cascading rules are not suitable for the reduced matrices [M'].There are two possible resolutions to this problem.First, one can develop a set of new cascading rules which preserves flux, cascade the matrices [M'] using these new rules, and compute internal fluxes for the devices.Alternatively, _one can transform the matrix [Iw'] to a Markov matrix [-&f], cascade the system using standard cascading rules (i.e., Ref. 7)) and once the fluxes are computed relate the internal fluxes computed from the [iM] to those fluxes that would have been computed by cascading ma&ix [-&$'I.
FIG. 2. Velocity field curve computed from the spectral flux method with three and five coetlicients.Also shown are the Monte Carlo results for FIG. 4. (a) The electric-field profile for the sample cakulation.(b) Average velocity xx position within the device.The solid line represents the Monte Carlo simulation results and the dotted line is from the simuiation using the spectral method.