A fast algorithm for the study of wave-packet scattering at disordered interfaces

We describe a very efficient algorithm for determining the scattering of wave-packets in 2D and 3D off a rough interface. The immediate objective is to derive effective scattering rates for utilisation in Monte Carlo simulation of MOSFET devices. Our methodology involves solving the multi-dimensional time-dependent Schrodinger equation, based on two-time iteration and direct integration.


INTRODUCTION
The study of the crucial interface roughness scattering in SiGe MOSFET devices is hampered by our lack of detailed models, particularly for situations where germanium desegregation and clustering occurs. Figure sketches a worse case scenario for a rough interface in which intercalation occurs so that the surface of constant potential is folded as compared to a simple random landscape model. Complex interface roughness scattering is common to all decanano silicon devices when the atomistic fluctuation potential due to discrete impurity/dopant distributions is encountered. Currently deployed models are somewhat simplistic and do not account for scattering in the atomistic domain [1,2]. The problem is essentially 3-D dimensional and it is computationally intensive for any modelling scheme. In the present paper we describe a very efficient algorithm for determining the scattering of wave-packets in 2D and 3D off a rough interface. The objective is to derive effective scattering rates for utilisation in Monte Carlo simulation of MOSFET devices. the encounter with the potential describing interface roughness will lead to scattering between propagating incoming states ks and outgoing states k at a transition rate R(k, k). This rate describes scattering off the entire interface and formally includes multiple scattering effects. Further transition rates may be defined which describe the temporary/permanent capture into the interface or through the interface by tunnelling and the reverse generation processes. To derive the transition rate consider an incident carrier in a plane wave state with wave vector ks which evolves into a state ](t;k)) describing the reflection/scattering of the carrier away from the interface. The incident wave is represented as the limit (L oo) of a normalised wave packet I(t;k,L)) of spatial scale L with mean wave vector ks. The probability amplitude that the scattered/reflected carrier makes a transition to a plane wave state k is the projection (k[(t,k)). It follows that the transition rate R(k, k) is the asymptotic time derivative in the limit (Lo): dC(ka,k;t) dl(kl(t,k.))l 2 R(k..k) dt dt (1) Practically, the transition rate may be calculated numerically by following the evolution of a wave packet [(t;k,,L)) incident on the interface and tracking the momentum probability distribution (2) Because the incident wave packet is finite and normalised the momentum distribution at any outgoing state k will be initially zero and then increase monotonically to a steady state P(k, o). The time derivative (d/dt)P(k, t) will correspondingly grow from zero, passing through a maximum (d/dt)P(k,t)ltmax corresponding to the escape of the centre of the wave packet; and finally it will return to zero, as the scattered wave packet becomes asymptotically free. In the case of a wide Gaussian packet, this maximum rate corresponds to the scattering rate for an incident plane wave. The scattering rate is finally identified as: In this derivation it is easy to formally calculate the scattering rate in terms of the T-matrix This fine-grained rate is usually averaged over a small neighbourhood of the exit wave vector to generate the coarse-grained scattering rates usually deployed in Monte Carlo. These depend explicitly on the energy density of final states p(sk): --- (7) in the limit of small At time interval. The Euler approximation is non-central in time, which rapidly leads to instabilities (loss of unitary) thus rendering it unsuitable for numerical simulations. The Cayley approximation used in the Crank-Nicholson method exhibits good numerical stability, but it is computationally expensive especially in higher-dimensions. For example, a 2-dimensional system, discretized into N spatial steps would require the inversion of an N 2 N 2 matrix. In our scheme we consider a two-time step iteration process at and t-+-At leading to: + t) (8) and hence we obtain the following finite-difference equation for at the next time step: In this scheme we must use two initial conditions one of which is prescribed as the initial state and the state at a later time step is determined by the Euler or Crank-Nicolson method. The complexity is only 30N 3 for 3D since this is a direct scheme and it requires no matrix inversion.

STABILITY ANALYSIS
The stability of the scheme may be analysed by a local yon Neumann procedure [6,7] where the independent solutions, Eigen modes of the difference equations are of the form: 2/) 7 (k) n e ikjAx (10) where k is a real spatial wave-number and (k) is complex number, referred to as the amplification factor for a given k. We have considered just onedimension with a constant effective mass, as these considerations are sufficient to derive the stability criterion that hold in general. Note that the largest wave-number that can be represented on a grid is k (1/Ax). For stability we require I(k)] _< 1.
Here the stability analysis gives a quadratic for (k), because of the occurrence of three-time steps and / At; in the finite difference approximation (Eq. (9)). Substitution for an eigenmode Eq. (10) into Eq. (9), gives: where Vmax is the maximum potential arising in the simulation. This condition is similar to the Courant condition, commonly arising in the numerical study of classical flux-conservative equations. Although the two-time method is not unconditionally stable like the Crank-Nicholson, this is not a limitation since the stability criterion are easily satisfied by the requirements (1/Ax) be greater than the maximum wave-vector we wish to resolve, and At is less than the smallest time scale that is to be considered. Indeed we would expect our method to exhibit better stability than the Crank-Nicholson, as less operations are required per time step reducing the possibility of rounding errors.

APPLICATION TO A ROUGH INTERFACE
The procedure may be illustrated by our preliminary applications to 2D surface roughness scattering on a fractal generated hard walled interface, which was generated using a random Weierstrass function: f(x) Z CkA-k sin (Akx + Ak) (14) k=l where Ck and Ak are independent random variables; Ck follows a Gaussian distribution, with mean,/z-0, and cr 2-1, whereas Ak is uniformly distributed on [0, 27r]. A (> 1) and a are chosen so as to insure that f(x) posses a r.m.s, vertical height deviation of A-2A, and an auto-correlation length of L 10 in agreement with published experimental data [8]. A typical slice through a rough surface is shown in Figure 2. Figure 3a shows the incident probability distribution (t-0) for position which gives rise to the distributions at a later time after scattering, shown for scattering off a smooth and rough interface in Figures 3b and 3c respectively. We note that Figure 3c shows strong streaming effects, this is because the rough interface distorts the momentum distribution parallel to the interface, while the smooth interface leaves this unchanged.  Smooth Interface FIGURE 3a An incident wave packet colliding with an interface on the right hand side of the picture, the incoming angle to the normal is r/4.

Smooth Interface
FIGURE 3b Outgoing wave packet having colliding with an interface on the right hand side of the picture, notice that the Gaussian shape is maintained due to specular reflection.
Rough Interface FIGURE 3c Outgoing wave packet having colliding with a rough interface as given by the random Weiestrass function. There is strong streaming in other directions other than just that due to specular reflection; a clear indication of interface roughness scattering. rates are appropriately modified and the free flight simulation then includes specular reflection at the walls. For soft walled model, or where capture or tunnelling is included this separation is not appropriate. In conclusion we have presented a new method for the evaluation of Monte Carlo scattering rates for complex 2-D and 3-D interface roughness models based on analysis of wave packet scattering. The direct integration algorithm presented is highly efficient for 2-D and 3-D and the stability is generally superior to Crank-Nicolson. 3 FIGURE 5 dP(k,t)/dt for a wave packet colliding off a rough interface. Notice that this rises to a maximum, which corresponds to the scattering rate generally used and while finally tends to zero, as the scattered wave packet escapes asymptotically free of the influence of the interface.

DISCUSSION AND CONCLUSIONS
For Monte Carlo applications it may useful to separate out the scattering from the spatially averaged hard interface, so that the effective scattering potential is V(r)-(V(r)). The exact