Solid-State Eleerronics zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA Vol. 32, No. 10, pp. 893-896, 1989 0038-l 101/89 $3.00 + 0.00 Printed in Great Britain. All rights reserved Copyright 0 1989 Pergamon Press plc zyxwvutsrq ITERATION APPROACH FOR SOLVING THE BOLTZMANN EQUATION WITH THE MONTE CARLO METHOD M. NEDJALKOV and P. VITANOV Institute of Microelectronics, Sofia, Bulgaria (Received 24 January 1989; in revised zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPO form 26 April 1989) Abstract-This study shows that the Monte Carlo technique of Jacoboni et al. can be obtained from the integral representation of the Boltzmann Transport Equation as well. A new iteration approach has been developed which allows problems connected with the enhancement of the computational efficiency, to be tackled. 1. INTRODUCTION The solution of the Boltzmann equation (BE) by the Monte Carlo (MC) method represents an established approach to semiconductor device modelling. In many cases, MC device models are generally accepted as containing a more rigorous description of device physics than deterministic models. The wide use of the MC method in semiconductor device simulation requires MC techniques optimized in achieving high accuracy together with low waste of processor time. The pioneer’s work appeared in 1985, when Socha and Krumhansl[l] created the path integral-an MC technique which guarantees high accuracy in the domain under consideration of the phase space. This technique requires knowledge of the distribution function in all moments previous to the current moment. In 1988 Jacoboni, Poli and Rota[Z] created the technique MC backwards. It enjoys the advantages of the path-integral, but here only knowledge is needed of the initial distribution function (for example in t = 0). It is established in Ref.[3] that the One particle, the Ensemble and the path-integral techniques though they appear quite different, have the same basethe integral representation of the Boltzmann equation (BE). In this paper it is shown that the technique of Jacoboni et al. can be obtained from the integral representation of the BE as well. The iteration ap- proach used allows: -consideration of the problems connected with the minimization of the variance, the possible error and the labor consumption-parameters controlling the computer waste time. -calculation of these parameters during the time of the distribution function simulation. 2. THE TECHNIQUE OF JACORONI, POLA AND ROTA This section considers the basic ideas developed by Jacoboni et al. The starting point is the standard form of BE (for simplicity a homogeneous system is con- sidered): where Si(k, k’)dk’ and S,(k,k’) dk’ are the same scattering frequencies for jumping from state k to states dk’ about zyxwvutsrqponmlkjihgfedcbaZYXWVU k’. The indices i-in and o-out clarify the balance of the particles incoming and outgoing from the classical trajectory. The time de- pendence of the wave vector is: k(f) = k - I , F(y ) dy so that k(r) = k. I’ After integration in the limits CW: f[k(t), tl= ‘dr, dk,Si(k,, k(t,If(k,, t,) s s zyxwvutsrqponmlkjihgfedcbaZYXWV 0 - 'dl, dk,W(r,h k,lf[k(t,), 41 s s 0 +fK’% 01. (2) The solution of this equation is found iteratively as a series f(k, t) =f(O) +f(‘) +f(2) +fu) + . After replacing f with zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQ k-l ,CoP (3) in (2) and grouping the terms with equal orders, the expression for f(‘) as a linear combination of the following integrals is obtained: I ss I dt, dkt$[kt,k(l,)l ss dt2 dk2 0 0 x %F,(t,)*k,l. . ..fW,(O), 01 (4) where the number of the multipliers in the combina- 893