Spherical Harmonic Modeling of a 0 . 05 Base BJT : A Comparison with Monte Carlo and Asymptotic Analysis

We perform a rigorous comparison between the Spherical Harmonic (SH) and Monte Carlo (MC) methods of solving the Boltzmann Transport Equation (BTE), on a 0.05 μm base BJT. We find the SH and the MC methods give very similar results for the energy distribution function, using an analytical band-structure, at all points within the tested devices. However, the SH method can be as much as seven thousand times faster than the MC approach for solving an identical problem. We explain the agreement by asymptotic analysis of the system of equations generated by the SH expansion of the BTE.


INTRODUCTION
The Spherical Harmonic (SH) method for solving the BTE has given rise to considerable contro- versy.Researchers developing the method have reported promising results [1][2][3].Other research-  ers have questioned whether the SH method can indeed provide the distribution function [4].Still  others have called for a study of various methods, including the SH, especially under conditions of rapid spatial variations [5].
In this paper, we respond to these questions by performing a rigorous comparison between the SH and MC methods of solving the Boltzmann equation, while using a 0.05 lm base BJT as a test vehicle.I,Ve find the SH and the MC methods give very similar results for the energy distribution function, using an analytical band-structure, at all points within the tested devices.However, the SH method can be as much as seven thousand times faster than the MC approach for solving an identical problem.
We explain the agreement between SH and MC by analysis of the SH system of equations.Asymptotic analyses show that, due to high scattering rates, high order SH terms decay quickly enough to have minimal effect on the energy distribution function, even for abrupt large spatial variations.Additionally, boundary condi- tions further restrict the population of higher Corresponding author.order terms, thereby facilitating use of a low order SH expansion.

OVERVIEW OF SPHERICAL HARMONIC METHOD
Since the main focus of this paper is to demon- strate and explain why the SH method can give relatively good agreement with MC simulations, we begin with an overview of the SH approach.First, we write the BTE below: (1) where (f') is the potential, f(', k) is the electron distribution function; k is the electron wave- vector; ' is the position vector; h is Planck's constant, and the sum is over the various collision processes represented by the transition rates Si.Now, we express the distribution function in terms of the following SH expansion: /=0 rn=-l here Y'(O,) are the spherical harmonic basis functions, 0 and are the polar and azimuthal angles of the wave-vector, respectively, fn(,,e) represents the expansion coefficients.The eventual goal of the approach is to determine the coeffi- cients and thereby the distribution function for the device.We showed previously that by substituting the spherical harmonics expansion for the dis- tribution function into the BTE, while taking advantage of the recurrence and orthogonality relationships between spherical harmonics the LHS of the BTE can be transformed into an infinite set of identical expressions for the coeffi- cients.We also showed the if we take scattering to be either isotropic or elastic (which is a reasonable approximation in silicon), and make use of the addition theorem, the collision integral can be transformed into a relatively simple expression [1, 2].Putting the LHS and the collision integral together gives the following equation for the zero order coefficient f0.0 1-7'() While the expressions for the coefficients for l> 0 are ,.. qEi(f') 0 2 i=1 0 (0 l+2-'/'()) where v(e) 2mf(e)/mf'(e), 7'(e) df(e)/de, 7(e) represents the dispersion relation, and Ei(?) is the electric field in the direction.The sum over represents the Cartesian directions of a 2-D device cross-section in real-space.The raising and low- ering operators fi{ which we developed to relate the coefficients are [1]: l,i(', )h(')d' In (5) we have introduced St,i, which is /'th coefficient of the Legendre polynomial expansion of the th type of scattering mechanism [2], as well as the density of states h(e) 4x/Tr/h3(m*3/271/2 (e)7'()).The superscript out refers to a process where the electron is scattered out of the state k to another state k', while the superscript in represents scattering into k.The superscripts q: in (5) refer to absorption (-) and emission (+) of energy respectively.(For elastic processes q: have been eliminated, since there is no energy emission or absorption.In cases where there is no process for energy emission (absorption), the superscript + (-) has been eliminated).
While a system of equations has been obtained to arbitrarily high order, we must truncate the system to actually obtain a solution.Since the scattering rate in silicon is relatively high, and silicon is a covalent semiconductor, a relatively low order truncation should be applicable.(This truncation will be verified later in the paper).We therefore use an expansion which includes the first four spherical harmonics (Y, Yi -1, Y, y0).After truncating, the total energy H is used as the independent energy variable as opposed to the kinetic energy e, (H=e-q) [3].We can then arrive at the followingtractable expression for the symmetrical part of the distribution function throughout the device [1, 2]: @/2(H + qdp)7-(H + qdp) VFOo(x H) -V 7y-i7(_ y, 3(o 1)On,, (H + q0)7'(H + 2h 3pn e,/KoT { + q, + q, n-+ #7(n + q* + x 7'(H+qO+hn)[en/KTF(x,y,H+hn) F (x, y, n)]) In (6) Dn optical phonon deformation potential; n optical phonon vibrational frequency; F0 symmetrical part of electron distribution; 1/-is the total phonon scattering rate.

BJT SIMULATION RESULTS AND DISCUSSION OF TRUNCATION
The main point of this paper is to demonstrate that a low order truncation is indeed suitable for determining the energy distribution function in semiclassical silicon devices under the spherical band approximation.To show this, we used the above formulation to simulate a 0.05 gm base BJT.
We chose the BJT because it represented a very aggressive challenge, with electric fields pointing in both directions, a narrow base of 50 nm, with fields varying as rapidly as 200 kV/cm over distances as small as 20 nm.As a benchmark, we also simulated the same device using the Monte Carlo approach.The input parameters (analytical band structure, intervalley, optical and acoustic phonon scattering) were identical for both the SH and MC calculations.Since the main objective was to determine if the low-order SH method could respond to the rapid field variations of the BJTs, we performed both the Monte Carlo and the Spherical Harmonic calculations as post-processor simulations.
We show some example results below.In Figures a and b we show the doping profile and resulting electric field with an applied poten- tial of VBE 1.0 V and VcB 3.0V.In Figures 2a   and 2b we show the 3-dimensional result for the energy distribution function calculated along the device.From a qualitative point of view, both distributions look alike, although the one gener- ated by the SH approach is populated over the complete range of the calculation whereas the MC results does not have data for high energies due to the time consuming stochastic nature of the MC process.In Figures 3a and 3b we show the energy distribution in more detail by plotting its values at regular intervals along the device.In the plots, we compare the SH results (solid lines)

FIGURE
The doping concentration (top) and electric field (bottom) for the 0.05 pm p-base Gaussian doped BJT with VBE 1.0 V and VcB 3.0 V.
. .results are in very good agreement for the entire distribution function.We find this to be an extremely encouraging result, considering that the SH method required only 10 sec to evaluate, while the MC approach took more than 10 hours on a DEC Alpha 266 MHz workstation.The large difference in CPU time is attributable to the deterministic nature of the SH approach, while the MC method requires considerable time to accumulate reasonable statistics for the high energy tail and to surmount the large potential barrier at the emitter-base junction.(The MC could probably be optimized but our attempts at statistical enhancement proved unreliable when subjected to scrutiny).
We explain the agreement between SH and MC by asymptotic analysis.First of all, for very large values of SH index l, l+l l, and Ft+l Under these conditions, Eq. ( 4) reduce to the following relatively simple expression.
v(H + qqb) OFn (x, H) Ft(x, H) ( Ox a(I-I + q) where A(H + qq) v(H+ qO)(H+ q(p) is a mean free path, the value of which depends on energy.
For large energy A becomes approximately con- stant, while for lower energies A can be approxi- mated as the following piecewise linear expression: A(H+ qb) c (H+ q) + /3 within [x0, x] for a constant H.By integrating (7) over position space from Xo to x for a constant H, using the chain rule dx-dxde, as well as approximating E(x) as constant in the small interval, we can obtain Fz(x, I F(xo, I-I) A(H + qbo) (8) which simplifies to the following expression in the high energy region where A is constant.
Ft(x, H) Ft(xo, H) exp ( x x0 ] (9) A ent in (8 becomes very large, and the term in the square brackets is always less than unity).For low to moderate energies, where there is a transition from a low electric field to a high one, or visa- versa, the monotonic nature of A again insures that the square bracketed term in (8) will be small.In addition to the asymptotic analysis for l, physical boundary conditions also suppress the magnitude of high order terms.For example, the physical requirement that the distribution function be single valued at the energy c 0 (is not dependent on angle at zero energy), leads to the requirement that all SH components with l> 0 must vanish at zero energy.Furthermore, the requirement that the distribution function is close to Maxwellian at the device contacts implies that higher order SH terms are negligible at these boundaries.Also, since in silicon the scattering rate is fairly high, and since it is a covalent semiconductor, we usually take the major scattering rates to be isotropic.This leads to scattering which tends to give rise to a distribution function which is largely spherical in nature.Finally, as one examines the generalized SH system given by Eqs. ( 3) and (4), it is clear that the equations are only coupled to nearest neighbors in l.It is therefore difficult to propagate information from low order terms to higher order ones, thereby minimizing the relative importance of higher order SH coefficients.
In high energy, (9) shows that F will decay exponentially in x-space with a characteristic length of A. Depending on the transport model one is using, A ranges from approximately 2 nm to 5 nm for energies greater than 1.5 eV in silicon.This indicates that high-energy spherical harmo- nics decay very quickly in comparison with typical dimensions in deep submicron devices.For regions of low energy and low electric field, (8) ensures that Ft also decays very rapidly with distance.
(While this is not immediately obvious, it is due to the relationship between 4 and E, and the fact that A is monotonically decreasing, ensures the expon-

FIGURE 2 FIGURE 3
FIGURE 2 The calculated electron distribution function by SH-BTE (top) and MC (bottom) for the 0.05 pm p-base Gaussian doped BJT.