- Email: [email protected]
Simulation of quantum transport by applying second-order differencing scheme to Wigner function model including spatially varying effective mass
PERGAMON
Solid-State Electronics 43 (1999) 81±86
Simulation of quantum transport by applying second-order dierencing scheme to Wigner function model including spatially varying eective mass Kyoung-Youm Kim, Byoungho Lee School of Electrical Engineering, Seoul National University, Kwanak-Gu, Shinlim-Dong, Seoul 151-742, South Korea Received 9 October 1997; received in revised form 7 May 1998
Abstract The discretization method and numerical calculations of the Wigner equation to introduce the in¯uence of spatially varying eective mass as well as to reduce the errors involved in the discretization are presented. The numerical errors due to a gradient operator for spatial derivative are considerable in the conventional upwind± downwind dierencing scheme. The errors can be reduced by adopting the second-order dierencing scheme. The eect of spatially varying eective mass is also very important and only by incorporating it simultaneously with the second-order dierencing scheme can we analyze the quantum devices exactly. # 1998 Elsevier Science Ltd. All rights reserved.
1. Introduction The progress in crystal growth technology makes it more important to simulate quantum transport precisely in quantum-size structures and devices whose behavior is dominated by quantum eects. The Wigner function model [1] has been used and proved to be useful in investigating the behavior of quantum devices, especially for resonant tunneling diodes (RTDs) [2±6]. This model enables us to calculate both small signal a.c. response and large-signal transient response as well as steady-state behavior in the form of current± voltage (I±V) characteristic curves. The I±V curves calculated from this model show the expected negative dierential-conductance regions, but the peak-to-valley current ratios are always smaller than those obtained from the static tunneling theory and are often less than those observed experimentally [7]. Tsuchiya et al. improved it by considering the spatially varying eective mass in the formulation of quantum Liouville equation [3]. In the work they used the usual ®rstorder upwind±downwind dierencing scheme (UDS) for the gradient operator calculation. However, analyzing the resonant tunneling assuming the eective mass
of electrons is constant, Buot and Jensen had shown that numerical errors involved in the UDS are quite considerable [4, 5]. They used second-order dierencing scheme (SDS) in which the spatial derivative is approximated by using a second-order interpolating polynomial. The peak-to-valley ratio calculated by the SDS Wigner function model is in quite good agreement with that obtained by the tunneling theory [8] with homogeneous eective mass. In this paper we will apply the SDS approach to the model by Tsuchiya et al. [3], in which spatially varying eective mass is included. To the authors' knowledge this is the ®rst trial to include in heterostructure quantum device simulation both the SDS approach and spatially varying eective mass eect.
2. Discretization of the Wigner equation The exact non-local Wigner representation of onedimensional (in x-direction) electron transport considering the spatially varying eective mass is written as [3]
0038-1101/98/$ - see front matter # 1998 Elsevier Science Ltd. All rights reserved. PII: S 0 0 3 8 - 1 1 0 1 ( 9 8 ) 0 0 2 0 1 - 9
82
K.-Y. Kim, B. Lee / Solid-State Electronics 43 (1999) 81±86
1 @ f x, k dk 0 @ f x, k 0 h M1 x, k ÿ k 0 @t @x ÿ1 8p 1 dk 0 0 ÿ h k M2 x, k ÿ k 0 f x, k 0 ÿ1 4p 1 dk 0 h M3 x, k ÿ k 0 ÿ1 16p 2 @ f x, k 0 02 0 ÿ 4k f x, k @x2 1 dk 0 0 @ f x, k 0 ÿ h k M4 x, k ÿ k 0 4p @x ÿ1 1 1 dk 0 0 0 ÿ V x, k ÿ k f x, k , h ÿ1 2p
1
where f(x, k) is the Wigner function de®ned as the Fourier transform of density matrix [2] f x, k
1 1 dx exp ÿjkxr x x, x ÿ x , 2 2
2:0
given in Ref. [3] (the authors of Ref. [3] implicitly assumed @f/@x = 0 in the reservoirs because there is a second-order derivative term in Eq. (1)), which were also used implicitly in Refs. [4, 5]. When the position coordinate x (0R xR L, where L is the device length) is represented by Nx mesh points with uniform mesh spacing Dx [=L/(Nx ÿ1)], the relative coordinate x is represented by Nx mesh points with uniform mesh spacing 2Dx. Then, if we take the transforms with respect to the relative coordinate x as in Eqs. (2.0), (2.1), (2.2), (2.3), (2.4) and (3), the resulting functions are periodic in k with a period of p/Dx [2]. Therefore, the k dimension (ÿp/2DxRkR p/2Dx) is discretized into Nk mesh points with uniform mesh spacing Dk = p/NkDx, i.e., k can take the values as follows. k 2 f p=Dx j ÿ 1=2=Nk ÿ 1=2,
j 1, 2, . . . , Nk g: 5
and
1 @ 1 M1 x, k ÿ k 0 2 dx sin k ÿ k 0 x, @ x mÿ x, x 0 2:1 1 @ 1 M2 x, k ÿ k 0 2 dx cos k ÿ k 0 x, @ x m x, x 0 2:2 1 sin k ÿ k 0 x M3 x, k ÿ k 0 2 dx , mÿ x, x 0
2:3
1 cos k ÿ k 0 x M4 x, k ÿ k 0 2 dx , m x, x 0
2:4
1 1 V x, k ÿ k 0 2 dx v x x 2 0 1 ÿ v x ÿ x sin k ÿ k 0 x, 2 1 1 1 ÿ , mÿ x, x m x x=2 m x ÿ x=2 1 m x,
x
1 1 : m x x=2 m x ÿ x=2
If Dx is taken to be the lattice constant (say, a) of the material of which the device is made, the k range (ÿp/2aR kRp2a) corresponds to only half of the ®rst Brillouin zone (ÿp/aR kR p/a). However, it is suciently large to contain almost all states of the particle under consideration in our problems (the actual model used in the simulation will be described in Section 3). Previous application [3] of Eq. (1) to RTD made use of the UDS for spatial derivative, which is accurate to order O(Dx) and given by @ f x, k f x2Dx , k ÿ f x, k 2 , @x Dx
6:0
4:0
where the signs refer to the upwind±downwind directions. In this paper we make use of the more accurate SDS, in which the spatial derivative at the desired point is approximated by using a second-order interpolating polynomial; it is accurate to O(Dx 2) and is given by [4, 5]
4:1
@ f x, k @x
3
In the above equations, n(x) and m(x) denote the potential energy and position-dependent eective mass, respectively. Details can be found in Ref. [3]. Numerical solutions to Eq. (1) are obtained using the ®nite-dierence method with boundary conditions
3f x, k ÿ 4f x2Dx , k f x22Dx , k 3 : 2Dx 6:1
Then, Eqs. (1), (2.0), (2.1), (2.2), (2.3), (2.4)±(4.0) and (4.1) are transformed into the following discretized equations:
K.-Y. Kim, B. Lee / Solid-State Electronics 43 (1999) 81±86
83
X M1 x, k ÿ k 0 @f h @t 2Dx k0 X ÿ3f x, k 0 4f x Dx , k 0 ÿ f x 2Dx , k 0 , k 0 < 0 ÿ h k 0 M2 x, k ÿ k 0 f x, k 0 0 0 0 0 3f x, k ÿ 4f x ÿ Dx , k f x ÿ 2Dx , k , k >0 k0 X M3 x, k ÿ k 0 X k 0 M4 x, k ÿ k 0 h f x Dx , k 0 ÿ 2 1 2D2x k 02 f x, k 0 f x ÿ Dx , k 0 ÿ h 2 2Dx Dx k0 k0 ÿ3f x, k 0 4f x Dx , k 0 ÿ f x 2Dx , k 0 , k 0 < 0 1X ÿ V x, k ÿ k 0 f x, k 0 , h k 0 3f x, k 0 ÿ 4f x ÿ Dx , k 0 f x ÿ 2Dx , k 0 , k 0 > 0 8 3 4 1 > > < mÿ x, k 0 ÿ mÿ x ÿ D , k 0 mÿ x ÿ 2D , k 0 , 1 X x x 0 0 M1 x, k ÿ k sin k ÿ k x > 4Nk Dx x 3 4 1 > :ÿ ÿ , mÿ x, k 0 mÿ x Dx , k 0 mÿ x 2Dx , k 0
9 k 0< 0 > > = , > ; k0 > 0 >
8 3 4 1 > ÿ , > < x, k 0 0 X m m m x ÿ D , k x ÿ 2Dx , k 0 1 x cos k ÿ k 0 x M2 x, k ÿ k 0 > 2Nk Dx x 3 4 1 > :ÿ ÿ , m x, k 0 m x Dx , k 0 m x 2Dx , k 0
9 k 0< 0 > > = , > ; k0 > 0 >
M3 x, k ÿ k 0
1 X sin k ÿ k 0 x , 4Nk x mÿ x, x
8:2
M4 x, k ÿ k 0
1 X cos k ÿ k 0 x , Nk x m x, x
8:3
2 X 1 V x, k ÿ k v x x Nk x 2 1 ÿ v x ÿ x sin k ÿ k 0 x: 2
Tsuchiya et al. pointed out that the number Nx is exactly limited to Nk/2 from the requirement that the discretized equations should be identical with those of the conventional Wigner equation when the eective mass is constant [3]. But in our approach when Nx is Nk/2, the o-diagonal elements of M4 do not vanish. Then our discretized equations (Eqs. (7), (8.0), (8.1), (8.2), (8.3) and (9)) are not identical with those of constant eective mass. Therefore we separate M4 into d.c. and a.c. parts as follows: 1 X cos k ÿ k 0 x M4 x, k ÿ k Nk x m x, x 1 X 1 1 ÿ cos k ÿ k 0 x Nk x m x, x Mc 1 X 1 cos k ÿ k 0 x , 10 Nk x Mc 0
cos k ÿ k 0 x
xmax X
x
9
8:0
8:1
where 1/Mc is the non-zero constant oset value of 1/m + (x, x) which is identical to 2/mGaAs in our heterostructure (the heterostructure will be described in Section 3). For the d.c. part we limit Nx to Nk/2 for k = k 0 , but Nk for k$ k 0 . For the a.c. part we limit Nx to Nk/2. The d.c. part for k$ k 0 vanishes because xmax X
0
7
x
p nk ÿ nk 0 cos 2Dx nx Dx Nk
Nk X 2p cos nk ÿ nk 0 nx 0, Nk nx
11
where k = p/Dx[(nk ÿ 0.5)/Nk ÿ 0.5] and x = 2Dxnx are used (nk and nx are positive integers). In this way, we can formulate the discrete equations for the Wigner equation considering both SDS and spatially varying eective mass. When the current densities are de®ned at mid-points between mesh points, we expect the discretized continuity equation to have the form @ n x 1 j x Dx =2 ÿ j x ÿ Dx =2 : ÿ @t e Dx
12
If the carrier density is de®ned as n x
X Dk k
2p
f x, k,
13
then the current density operator must be obtained in such a way that it satis®es the discretized continuity
84
K.-Y. Kim, B. Lee / Solid-State Electronics 43 (1999) 81±86
equation (Eq. (12)) which is derived from the discretized Wigner equation (Eq. (7)) by summing over k. However, in the SDS approach, the eective mass contribution to the current operator is not local. The current operator should be given as follows in the nonlocal form: 1 eDk X 0 j x Dx h k M4 x, k ÿ k 0 2 4p k,k 0 3f x Dx , k 0 ÿ f x 2Dx , k 0 , k 0 < 0 : 3f x, k 0 ÿ f x ÿ Dx , k 0 , k0 > 0 14
3. Comparison between the SDS method and the UDS method The RTD model used in the simulation is shown in Fig. 1. It consists of a 4.5-nm-wide GaAs quantum well and 2.8-nm-wide Al0.3Ga0.7As barrier layers. The GaAs electrode layers of 17.5 nm are included on each side of the device. Outside the electrodes, ideal reservoirs characterized by the thermal equilibrium distributions of carriers are assumed. The conduction band discontinuity is taken to be 0.27 eV [2, 3]. Nx = 82, Nk = 60 and the mesh spacing Dx = 0.565 nm are used. The mesh spacing we use here is comparable to the lattice constant of GaAs [2]. The eective masses of electrons are 0.067m0 in GaAs and 0.092m0 in Al0.3Ga0.7As, respectively, where m0 is the free electron mass [3]. The doping density in the GaAs electrodes is given as 2 1018 cm ÿ 3 and the barriers and well layers are assumed to be undoped [2, 3]. Scattering processes are neglected and all calculations are performed at a room temperature of 300 K [2, 3]. I±V characteristics calculated by both the SDS and UDS with constant and spatially varying eective masses are shown in Fig. 2. The peak-to-valley ratio predicted by the tunneling theory [7] is about 9:1 (scat-
Fig. 1. Resonant tunneling diode model used in the simulation.
tering is ignored) and experimental data are above 5:1, especially for low temperature cases [9, 10] where phonon scattering does not play great roles. It is well known that phonon scattering process reduces the peak-to-valley ratio greatly [2, 4±7, 11, 12] and that when scattering is ignored the I±V curves at 300 and 77 K are almost the same, especially for the low bias regions where tunneling process dominates, but the peak-to-valley ratio of I±V curve at 300 K is nearly half of that at 77 K mainly due to the temperature-sensitive valley current (the valley current is the contribution of thermally excited electrons) [12]. Jensen and Buot show that when phonon scattering is considered, the peak-to-valley ratio at 77 K decreases by a factor of 2 [12], which suggests that the peak-to-valley ratio of the I±V curve at 300 K with scattering ignored is nearly the same as that at 77 K with scattering included (i.e., the experimental value). Therefore in our problem, although we perform simulation at room temperature, because we exclude scattering eect, the calculation must somehow predict a peak-to-valley ratio close to low temperature experimental values. Ray et al. [9] observed the peak-to-valley ratio of 6:1 at 77 K with the structure of 5-nm-thick GaAs layer sandwiched between two 5-nm-thick Al0.45Ga0.55As layers with the doping density at GaAs electrode of 2 1018 cm ÿ 3. Huang et al. [10] observed the peak-tovalley ratio of 6.2±7.7:1 at 77 K with the structure of 5-nm-thick GaAs layer sandwiched between two 5-nmthick Al0.3Ga0.7As layers with the doping density at GaAs electrode of 2 1018 cm ÿ 3. In both experiments, the peak-to-valley ratios are reduced greatly at 300 K
Fig. 2. Comparison between the I±V characteristics calculated by the UDS and SDS methods with constant and spatially varying eective masses. Starred (*) numbers show peak-tovalley ratios. (1) SDS calculation with spatially varying eective mass, (2) SDS calculation with constant eective mass m* = 0.067m0, (3) UDS calculation with spatially varying eective mass and (4) SDS calculation with constant eective mass m* = 0.067m0.
K.-Y. Kim, B. Lee / Solid-State Electronics 43 (1999) 81±86
85
(2.5:1 [9], 1.85±2.2:1 [10]) due to phonon scattering. The RTD model structure used in this paper is dierent from above experimental ones in barrier width. However, the barrier width change gives only a little variation in the peak-to-valley ratio [11]; it only takes much time and memory in the simulation. Therefore, our simulation model should be able to predict the peak-to-valley ratio close to the above experimental values at low temperature. Calculated peak-to-valley ratios are also shown in Fig. 2. As is apparent from Fig. 2, the spatially varying eective mass plays a great role in the accurate analysis of the device. But, the improvement of calculation by adopting SDS is greater than the eect of spatially varying eective mass. In the SDS calculation, if we consider the eect of spatially varying eective mass, the peak bias is not altered but the valley bias is decreased by 40 meV. The peak-to-valley ratio is 8.8:1, much larger than that of 4.4:1 obtained assuming the homogeneous eective mass, and very close to the experimental values. Therefore, under the Wigner transport formalism both the SDS approach and the inclusion of spatially varying eective mass are needed for more exact calculation of device characteristics. The reason why the spatially varying eective mass improves the peak-to-valley ratio can be explained as follows. In the tunneling current regions (low bias), the quantum correlation of potential term (Eq. (3)) is very important [2]. The tunneling phenomenon can be described by this non-local potential term. When the spatially varying eective mass is introduced, there is another kind of quantum correlation of eective mass (Eqs. (2.1), (2.2), (2.3) and (2.4), which are also very important in the low bias regions) and this correlation seems to enhance the tunneling phenomenon.
However, in the classical regions (high bias), the current density is proportional to the inverse of eective mass (see Eqs. (4.1), (8.3) and (14)). Therefore, the heavier (compared with homogeneous eective mass case) eective mass in the barrier makes the valley current density smaller. Due to the above two reasons, the peak-to-valley ratio is improved if we include the spatially varying eective mass eect. That is, the constant eective mass model underestimates the peak current and overestimates the currents in classical regions. It is possible to extend the simulation so that it includes the Coulomb interaction eect by self-consistent potential which can be obtained by solving Poisson's equation using the method of Ref. [3]. In the simulation, the iterative calculation stops when the change of potential energy becomes less than 10 meV at any position. Fig. 3 shows the calculated potential energy and electron density when the bias voltage is 0.26 V (peak point, dashed line) and 0.36 V (valley point, solid line), respectively. With the peak bias, we can see electron accumulation in the well and with the valley bias we can see electron accumulation in front of the ®rst potential barrier. For more accurate analysis, inelastic processes such as electron±electron collision or electron±phonon interaction must be included in the simulation [4, 5, 12]. Fig. 4 is the plot of I±V characteristics for both the ¯at-band model and self-consistent potential model. The scattering process is neglected. The peak and valley biases of the self-consistent potential model are higher than those of the ¯at-band model by about 100±130 meV. It is because the net bias forced on the quantum well becomes smaller when we consider the self-consistent potential due to the bias contributing to the bending of electrode potential [compare Fig. 3(a) with Fig. 1].
Fig. 3. Self-consistent potential (a) and electron density (b) in the RTD. Solid line is for the bias of valley voltage (0.36 V) and dashed line is for the peak voltage (0.26 V).
Fig. 4. Comparison between the I±V characteristics calculated by the ¯at-band model (dashed line) and the self-consistent potential model (solid line).
86
K.-Y. Kim, B. Lee / Solid-State Electronics 43 (1999) 81±86
4. Conclusion
References
The discretization method and numerical calculations of the Wigner equation are presented to include the eect of spatially varying eective mass as well as to reduce the errors involved in the discrete form of the equation. The numerical errors due to the gradient operator for spatial derivative are considerable in the conventional UDS method. The errors can be reduced by adopting the SDS method. The eect of spatially varying eective mass is very important in the quantum device simulation. Only by incorporating it simultaneously with the second-order dierencing scheme can we predict the experimental results more exactly.
[1] Wigner E. Phys Rev 1932;40:749. [2] Frensley WR. Phys Rev B 1987;36:1570. [3] Tsuchiya H, Ogawa M, Miyoshi T. IEEE Trans Electron Devices 1991;38:1246. [4] Jensen KL, Buot FA. J Appl Phys 1990;67:2153. [5] Buot FA, Jensen KL. Phys Rev B 1990;42:9429. [6] Frensley WR. Rev Mod Phys 1990;62:745. [7] Frensley WR. In: Einspruch NG, Frensley WR, editors. Heterostructures and Quantum Devices, New York: Academic Press, 1994. Chap. 9. [8] Hellman ES, Lear KL, Harris JS, Jr. J Appl Phys 1988;64:2798. [9] Ray S, Ruden P, Sokolov V, Kolbas R, Boonstra T, Williams J. Appl Phys Lett 1986;48:1666. [10] Huang CI, Paulus MJ, Bozada CA, Dudley SC, Evans KR, Stutz CE, Jones RL, Cheney ME. Appl Phys Lett 1987;51:121. [11] Tsuchiya H, Ogawa M, Miyoshi T. Jpn J Appl Phys 1992;31:745. [12] Jensen KL, Buot FA. IEEE Trans Electron Devices 1991;38:2337.
Acknowledgements The authors acknowledge the support of the Ministry of Science and Technology of Korea through Basic Future Technology Project of MOST.
Solid-State Electronics 43 (1999) 81±86
Simulation of quantum transport by applying second-order dierencing scheme to Wigner function model including spatially varying eective mass Kyoung-Youm Kim, Byoungho Lee School of Electrical Engineering, Seoul National University, Kwanak-Gu, Shinlim-Dong, Seoul 151-742, South Korea Received 9 October 1997; received in revised form 7 May 1998
Abstract The discretization method and numerical calculations of the Wigner equation to introduce the in¯uence of spatially varying eective mass as well as to reduce the errors involved in the discretization are presented. The numerical errors due to a gradient operator for spatial derivative are considerable in the conventional upwind± downwind dierencing scheme. The errors can be reduced by adopting the second-order dierencing scheme. The eect of spatially varying eective mass is also very important and only by incorporating it simultaneously with the second-order dierencing scheme can we analyze the quantum devices exactly. # 1998 Elsevier Science Ltd. All rights reserved.
1. Introduction The progress in crystal growth technology makes it more important to simulate quantum transport precisely in quantum-size structures and devices whose behavior is dominated by quantum eects. The Wigner function model [1] has been used and proved to be useful in investigating the behavior of quantum devices, especially for resonant tunneling diodes (RTDs) [2±6]. This model enables us to calculate both small signal a.c. response and large-signal transient response as well as steady-state behavior in the form of current± voltage (I±V) characteristic curves. The I±V curves calculated from this model show the expected negative dierential-conductance regions, but the peak-to-valley current ratios are always smaller than those obtained from the static tunneling theory and are often less than those observed experimentally [7]. Tsuchiya et al. improved it by considering the spatially varying eective mass in the formulation of quantum Liouville equation [3]. In the work they used the usual ®rstorder upwind±downwind dierencing scheme (UDS) for the gradient operator calculation. However, analyzing the resonant tunneling assuming the eective mass
of electrons is constant, Buot and Jensen had shown that numerical errors involved in the UDS are quite considerable [4, 5]. They used second-order dierencing scheme (SDS) in which the spatial derivative is approximated by using a second-order interpolating polynomial. The peak-to-valley ratio calculated by the SDS Wigner function model is in quite good agreement with that obtained by the tunneling theory [8] with homogeneous eective mass. In this paper we will apply the SDS approach to the model by Tsuchiya et al. [3], in which spatially varying eective mass is included. To the authors' knowledge this is the ®rst trial to include in heterostructure quantum device simulation both the SDS approach and spatially varying eective mass eect.
2. Discretization of the Wigner equation The exact non-local Wigner representation of onedimensional (in x-direction) electron transport considering the spatially varying eective mass is written as [3]
0038-1101/98/$ - see front matter # 1998 Elsevier Science Ltd. All rights reserved. PII: S 0 0 3 8 - 1 1 0 1 ( 9 8 ) 0 0 2 0 1 - 9
82
K.-Y. Kim, B. Lee / Solid-State Electronics 43 (1999) 81±86
1 @ f x, k dk 0 @ f x, k 0 h M1 x, k ÿ k 0 @t @x ÿ1 8p 1 dk 0 0 ÿ h k M2 x, k ÿ k 0 f x, k 0 ÿ1 4p 1 dk 0 h M3 x, k ÿ k 0 ÿ1 16p 2 @ f x, k 0 02 0 ÿ 4k f x, k @x2 1 dk 0 0 @ f x, k 0 ÿ h k M4 x, k ÿ k 0 4p @x ÿ1 1 1 dk 0 0 0 ÿ V x, k ÿ k f x, k , h ÿ1 2p
1
where f(x, k) is the Wigner function de®ned as the Fourier transform of density matrix [2] f x, k
1 1 dx exp ÿjkxr x x, x ÿ x , 2 2
2:0
given in Ref. [3] (the authors of Ref. [3] implicitly assumed @f/@x = 0 in the reservoirs because there is a second-order derivative term in Eq. (1)), which were also used implicitly in Refs. [4, 5]. When the position coordinate x (0R xR L, where L is the device length) is represented by Nx mesh points with uniform mesh spacing Dx [=L/(Nx ÿ1)], the relative coordinate x is represented by Nx mesh points with uniform mesh spacing 2Dx. Then, if we take the transforms with respect to the relative coordinate x as in Eqs. (2.0), (2.1), (2.2), (2.3), (2.4) and (3), the resulting functions are periodic in k with a period of p/Dx [2]. Therefore, the k dimension (ÿp/2DxRkR p/2Dx) is discretized into Nk mesh points with uniform mesh spacing Dk = p/NkDx, i.e., k can take the values as follows. k 2 f p=Dx j ÿ 1=2=Nk ÿ 1=2,
j 1, 2, . . . , Nk g: 5
and
1 @ 1 M1 x, k ÿ k 0 2 dx sin k ÿ k 0 x, @ x mÿ x, x 0 2:1 1 @ 1 M2 x, k ÿ k 0 2 dx cos k ÿ k 0 x, @ x m x, x 0 2:2 1 sin k ÿ k 0 x M3 x, k ÿ k 0 2 dx , mÿ x, x 0
2:3
1 cos k ÿ k 0 x M4 x, k ÿ k 0 2 dx , m x, x 0
2:4
1 1 V x, k ÿ k 0 2 dx v x x 2 0 1 ÿ v x ÿ x sin k ÿ k 0 x, 2 1 1 1 ÿ , mÿ x, x m x x=2 m x ÿ x=2 1 m x,
x
1 1 : m x x=2 m x ÿ x=2
If Dx is taken to be the lattice constant (say, a) of the material of which the device is made, the k range (ÿp/2aR kRp2a) corresponds to only half of the ®rst Brillouin zone (ÿp/aR kR p/a). However, it is suciently large to contain almost all states of the particle under consideration in our problems (the actual model used in the simulation will be described in Section 3). Previous application [3] of Eq. (1) to RTD made use of the UDS for spatial derivative, which is accurate to order O(Dx) and given by @ f x, k f x2Dx , k ÿ f x, k 2 , @x Dx
6:0
4:0
where the signs refer to the upwind±downwind directions. In this paper we make use of the more accurate SDS, in which the spatial derivative at the desired point is approximated by using a second-order interpolating polynomial; it is accurate to O(Dx 2) and is given by [4, 5]
4:1
@ f x, k @x
3
In the above equations, n(x) and m(x) denote the potential energy and position-dependent eective mass, respectively. Details can be found in Ref. [3]. Numerical solutions to Eq. (1) are obtained using the ®nite-dierence method with boundary conditions
3f x, k ÿ 4f x2Dx , k f x22Dx , k 3 : 2Dx 6:1
Then, Eqs. (1), (2.0), (2.1), (2.2), (2.3), (2.4)±(4.0) and (4.1) are transformed into the following discretized equations:
K.-Y. Kim, B. Lee / Solid-State Electronics 43 (1999) 81±86
83
X M1 x, k ÿ k 0 @f h @t 2Dx k0 X ÿ3f x, k 0 4f x Dx , k 0 ÿ f x 2Dx , k 0 , k 0 < 0 ÿ h k 0 M2 x, k ÿ k 0 f x, k 0 0 0 0 0 3f x, k ÿ 4f x ÿ Dx , k f x ÿ 2Dx , k , k >0 k0 X M3 x, k ÿ k 0 X k 0 M4 x, k ÿ k 0 h f x Dx , k 0 ÿ 2 1 2D2x k 02 f x, k 0 f x ÿ Dx , k 0 ÿ h 2 2Dx Dx k0 k0 ÿ3f x, k 0 4f x Dx , k 0 ÿ f x 2Dx , k 0 , k 0 < 0 1X ÿ V x, k ÿ k 0 f x, k 0 , h k 0 3f x, k 0 ÿ 4f x ÿ Dx , k 0 f x ÿ 2Dx , k 0 , k 0 > 0 8 3 4 1 > > < mÿ x, k 0 ÿ mÿ x ÿ D , k 0 mÿ x ÿ 2D , k 0 , 1 X x x 0 0 M1 x, k ÿ k sin k ÿ k x > 4Nk Dx x 3 4 1 > :ÿ ÿ , mÿ x, k 0 mÿ x Dx , k 0 mÿ x 2Dx , k 0
9 k 0< 0 > > = , > ; k0 > 0 >
8 3 4 1 > ÿ , > < x, k 0 0 X m m m x ÿ D , k x ÿ 2Dx , k 0 1 x cos k ÿ k 0 x M2 x, k ÿ k 0 > 2Nk Dx x 3 4 1 > :ÿ ÿ , m x, k 0 m x Dx , k 0 m x 2Dx , k 0
9 k 0< 0 > > = , > ; k0 > 0 >
M3 x, k ÿ k 0
1 X sin k ÿ k 0 x , 4Nk x mÿ x, x
8:2
M4 x, k ÿ k 0
1 X cos k ÿ k 0 x , Nk x m x, x
8:3
2 X 1 V x, k ÿ k v x x Nk x 2 1 ÿ v x ÿ x sin k ÿ k 0 x: 2
Tsuchiya et al. pointed out that the number Nx is exactly limited to Nk/2 from the requirement that the discretized equations should be identical with those of the conventional Wigner equation when the eective mass is constant [3]. But in our approach when Nx is Nk/2, the o-diagonal elements of M4 do not vanish. Then our discretized equations (Eqs. (7), (8.0), (8.1), (8.2), (8.3) and (9)) are not identical with those of constant eective mass. Therefore we separate M4 into d.c. and a.c. parts as follows: 1 X cos k ÿ k 0 x M4 x, k ÿ k Nk x m x, x 1 X 1 1 ÿ cos k ÿ k 0 x Nk x m x, x Mc 1 X 1 cos k ÿ k 0 x , 10 Nk x Mc 0
cos k ÿ k 0 x
xmax X
x
9
8:0
8:1
where 1/Mc is the non-zero constant oset value of 1/m + (x, x) which is identical to 2/mGaAs in our heterostructure (the heterostructure will be described in Section 3). For the d.c. part we limit Nx to Nk/2 for k = k 0 , but Nk for k$ k 0 . For the a.c. part we limit Nx to Nk/2. The d.c. part for k$ k 0 vanishes because xmax X
0
7
x
p nk ÿ nk 0 cos 2Dx nx Dx Nk
Nk X 2p cos nk ÿ nk 0 nx 0, Nk nx
11
where k = p/Dx[(nk ÿ 0.5)/Nk ÿ 0.5] and x = 2Dxnx are used (nk and nx are positive integers). In this way, we can formulate the discrete equations for the Wigner equation considering both SDS and spatially varying eective mass. When the current densities are de®ned at mid-points between mesh points, we expect the discretized continuity equation to have the form @ n x 1 j x Dx =2 ÿ j x ÿ Dx =2 : ÿ @t e Dx
12
If the carrier density is de®ned as n x
X Dk k
2p
f x, k,
13
then the current density operator must be obtained in such a way that it satis®es the discretized continuity
84
K.-Y. Kim, B. Lee / Solid-State Electronics 43 (1999) 81±86
equation (Eq. (12)) which is derived from the discretized Wigner equation (Eq. (7)) by summing over k. However, in the SDS approach, the eective mass contribution to the current operator is not local. The current operator should be given as follows in the nonlocal form: 1 eDk X 0 j x Dx h k M4 x, k ÿ k 0 2 4p k,k 0 3f x Dx , k 0 ÿ f x 2Dx , k 0 , k 0 < 0 : 3f x, k 0 ÿ f x ÿ Dx , k 0 , k0 > 0 14
3. Comparison between the SDS method and the UDS method The RTD model used in the simulation is shown in Fig. 1. It consists of a 4.5-nm-wide GaAs quantum well and 2.8-nm-wide Al0.3Ga0.7As barrier layers. The GaAs electrode layers of 17.5 nm are included on each side of the device. Outside the electrodes, ideal reservoirs characterized by the thermal equilibrium distributions of carriers are assumed. The conduction band discontinuity is taken to be 0.27 eV [2, 3]. Nx = 82, Nk = 60 and the mesh spacing Dx = 0.565 nm are used. The mesh spacing we use here is comparable to the lattice constant of GaAs [2]. The eective masses of electrons are 0.067m0 in GaAs and 0.092m0 in Al0.3Ga0.7As, respectively, where m0 is the free electron mass [3]. The doping density in the GaAs electrodes is given as 2 1018 cm ÿ 3 and the barriers and well layers are assumed to be undoped [2, 3]. Scattering processes are neglected and all calculations are performed at a room temperature of 300 K [2, 3]. I±V characteristics calculated by both the SDS and UDS with constant and spatially varying eective masses are shown in Fig. 2. The peak-to-valley ratio predicted by the tunneling theory [7] is about 9:1 (scat-
Fig. 1. Resonant tunneling diode model used in the simulation.
tering is ignored) and experimental data are above 5:1, especially for low temperature cases [9, 10] where phonon scattering does not play great roles. It is well known that phonon scattering process reduces the peak-to-valley ratio greatly [2, 4±7, 11, 12] and that when scattering is ignored the I±V curves at 300 and 77 K are almost the same, especially for the low bias regions where tunneling process dominates, but the peak-to-valley ratio of I±V curve at 300 K is nearly half of that at 77 K mainly due to the temperature-sensitive valley current (the valley current is the contribution of thermally excited electrons) [12]. Jensen and Buot show that when phonon scattering is considered, the peak-to-valley ratio at 77 K decreases by a factor of 2 [12], which suggests that the peak-to-valley ratio of the I±V curve at 300 K with scattering ignored is nearly the same as that at 77 K with scattering included (i.e., the experimental value). Therefore in our problem, although we perform simulation at room temperature, because we exclude scattering eect, the calculation must somehow predict a peak-to-valley ratio close to low temperature experimental values. Ray et al. [9] observed the peak-to-valley ratio of 6:1 at 77 K with the structure of 5-nm-thick GaAs layer sandwiched between two 5-nm-thick Al0.45Ga0.55As layers with the doping density at GaAs electrode of 2 1018 cm ÿ 3. Huang et al. [10] observed the peak-tovalley ratio of 6.2±7.7:1 at 77 K with the structure of 5-nm-thick GaAs layer sandwiched between two 5-nmthick Al0.3Ga0.7As layers with the doping density at GaAs electrode of 2 1018 cm ÿ 3. In both experiments, the peak-to-valley ratios are reduced greatly at 300 K
Fig. 2. Comparison between the I±V characteristics calculated by the UDS and SDS methods with constant and spatially varying eective masses. Starred (*) numbers show peak-tovalley ratios. (1) SDS calculation with spatially varying eective mass, (2) SDS calculation with constant eective mass m* = 0.067m0, (3) UDS calculation with spatially varying eective mass and (4) SDS calculation with constant eective mass m* = 0.067m0.
K.-Y. Kim, B. Lee / Solid-State Electronics 43 (1999) 81±86
85
(2.5:1 [9], 1.85±2.2:1 [10]) due to phonon scattering. The RTD model structure used in this paper is dierent from above experimental ones in barrier width. However, the barrier width change gives only a little variation in the peak-to-valley ratio [11]; it only takes much time and memory in the simulation. Therefore, our simulation model should be able to predict the peak-to-valley ratio close to the above experimental values at low temperature. Calculated peak-to-valley ratios are also shown in Fig. 2. As is apparent from Fig. 2, the spatially varying eective mass plays a great role in the accurate analysis of the device. But, the improvement of calculation by adopting SDS is greater than the eect of spatially varying eective mass. In the SDS calculation, if we consider the eect of spatially varying eective mass, the peak bias is not altered but the valley bias is decreased by 40 meV. The peak-to-valley ratio is 8.8:1, much larger than that of 4.4:1 obtained assuming the homogeneous eective mass, and very close to the experimental values. Therefore, under the Wigner transport formalism both the SDS approach and the inclusion of spatially varying eective mass are needed for more exact calculation of device characteristics. The reason why the spatially varying eective mass improves the peak-to-valley ratio can be explained as follows. In the tunneling current regions (low bias), the quantum correlation of potential term (Eq. (3)) is very important [2]. The tunneling phenomenon can be described by this non-local potential term. When the spatially varying eective mass is introduced, there is another kind of quantum correlation of eective mass (Eqs. (2.1), (2.2), (2.3) and (2.4), which are also very important in the low bias regions) and this correlation seems to enhance the tunneling phenomenon.
However, in the classical regions (high bias), the current density is proportional to the inverse of eective mass (see Eqs. (4.1), (8.3) and (14)). Therefore, the heavier (compared with homogeneous eective mass case) eective mass in the barrier makes the valley current density smaller. Due to the above two reasons, the peak-to-valley ratio is improved if we include the spatially varying eective mass eect. That is, the constant eective mass model underestimates the peak current and overestimates the currents in classical regions. It is possible to extend the simulation so that it includes the Coulomb interaction eect by self-consistent potential which can be obtained by solving Poisson's equation using the method of Ref. [3]. In the simulation, the iterative calculation stops when the change of potential energy becomes less than 10 meV at any position. Fig. 3 shows the calculated potential energy and electron density when the bias voltage is 0.26 V (peak point, dashed line) and 0.36 V (valley point, solid line), respectively. With the peak bias, we can see electron accumulation in the well and with the valley bias we can see electron accumulation in front of the ®rst potential barrier. For more accurate analysis, inelastic processes such as electron±electron collision or electron±phonon interaction must be included in the simulation [4, 5, 12]. Fig. 4 is the plot of I±V characteristics for both the ¯at-band model and self-consistent potential model. The scattering process is neglected. The peak and valley biases of the self-consistent potential model are higher than those of the ¯at-band model by about 100±130 meV. It is because the net bias forced on the quantum well becomes smaller when we consider the self-consistent potential due to the bias contributing to the bending of electrode potential [compare Fig. 3(a) with Fig. 1].
Fig. 3. Self-consistent potential (a) and electron density (b) in the RTD. Solid line is for the bias of valley voltage (0.36 V) and dashed line is for the peak voltage (0.26 V).
Fig. 4. Comparison between the I±V characteristics calculated by the ¯at-band model (dashed line) and the self-consistent potential model (solid line).
86
K.-Y. Kim, B. Lee / Solid-State Electronics 43 (1999) 81±86
4. Conclusion
References
The discretization method and numerical calculations of the Wigner equation are presented to include the eect of spatially varying eective mass as well as to reduce the errors involved in the discrete form of the equation. The numerical errors due to the gradient operator for spatial derivative are considerable in the conventional UDS method. The errors can be reduced by adopting the SDS method. The eect of spatially varying eective mass is very important in the quantum device simulation. Only by incorporating it simultaneously with the second-order dierencing scheme can we predict the experimental results more exactly.
[1] Wigner E. Phys Rev 1932;40:749. [2] Frensley WR. Phys Rev B 1987;36:1570. [3] Tsuchiya H, Ogawa M, Miyoshi T. IEEE Trans Electron Devices 1991;38:1246. [4] Jensen KL, Buot FA. J Appl Phys 1990;67:2153. [5] Buot FA, Jensen KL. Phys Rev B 1990;42:9429. [6] Frensley WR. Rev Mod Phys 1990;62:745. [7] Frensley WR. In: Einspruch NG, Frensley WR, editors. Heterostructures and Quantum Devices, New York: Academic Press, 1994. Chap. 9. [8] Hellman ES, Lear KL, Harris JS, Jr. J Appl Phys 1988;64:2798. [9] Ray S, Ruden P, Sokolov V, Kolbas R, Boonstra T, Williams J. Appl Phys Lett 1986;48:1666. [10] Huang CI, Paulus MJ, Bozada CA, Dudley SC, Evans KR, Stutz CE, Jones RL, Cheney ME. Appl Phys Lett 1987;51:121. [11] Tsuchiya H, Ogawa M, Miyoshi T. Jpn J Appl Phys 1992;31:745. [12] Jensen KL, Buot FA. IEEE Trans Electron Devices 1991;38:2337.
Acknowledgements The authors acknowledge the support of the Ministry of Science and Technology of Korea through Basic Future Technology Project of MOST.