Mô phỏng ba chiều linh kiện na-nô bán dẫn với lời giải phương trình poisson dựa trên thuật toán Gpbicg

215
TẠP CHÍ KHOA HỌC, Đại học Huế, Số 65, 2011 
MÔ PHỎNG BA CHIỀU LINH KIỆN NA-NÔ BÁN DẪN 
VỚI LỜI GIẢI PHƯƠNG TRÌNH POISSON 
DỰA TRÊN THUẬT TOÁN GPBICG 
Đinh Như Thảo, Dương Thị Diễm My, 
Nguyễn Châu Phương Thi, Ngô Thanh Thủy 
 Trường Đại học Sư phạm, Đại học Huế 
TÓM TẮT 
Bài báo trình bày việc xây dựng chương trình giải phương trình Poisson ba chiều dựa 
trên thuật toán GPBICG để sử dụng trong chương trình mô phỏng linh kiện na-nô bán dẫn bằng 
phương pháp Monte – Carlo tập hợp tự hợp. Chương trình mô phỏng được áp dụng để mô 
phỏng động lực học ba chiều của hạt tải trong các đi-ốt p-i-n bán dẫn GaAs. Các kết quả mô 
phỏng thu được hoàn toàn phù hợp với các kết quả của các công trình đã được công bố trước 
đây [1, 2]. Các kết quả chỉ ra rằng, chương trình giải phương trình Poisson dựa trên thuật toán 
GPBICG không những có tốc độ hội tụ nhanh mà còn có tính ổn định cao hơn các chương trình 
từng được sử dụng [2]. 
Từ khóa: Mô phỏng linh kiện bán dẫn, phương trình Poisson ba chiều, thuật toán 
GPBICG, phương pháp Monte – Carlo. 
1. Giới thiệu 
Nghiên cứu và phát triển các linh kiện na-nô bán dẫn đang thu hút sự quan tâm 
mạnh mẽ của giới khoa học do tính ứng dụng cao của nó [1, 2, 3]. Nghiên cứu thực 
nghiệm các linh kiện na-nô nói chung là rất tốn kém, đòi hỏi phải sử dụng công nghệ 
cao và mất nhiều thời gian. Các phương pháp nghiên cứu lý thuyết có thể giúp khắc 
phục được các hạn chế nêu ở trên, đặc biệt là các phương pháp mô phỏng như: phương 
pháp Monte – Carlo tập hợp tự hợp, phương pháp các phương trình cân bằng, mô hình 
kéo theo – khuếch tán [4]. Trong lớp các phương pháp đó, Monte – Carlo tập hợp tự hợp 
có nhiều ưu điểm nổi trội đặc biệt là tính chính xác và tính ổn định. Đây là phương pháp 
bán cổ điển với tốc độ tán xạ được tính toán dựa trên qui tắc vàng Fermi và việc khảo 
sát động lực học của hạt tải dựa trên các phương trình động học của Newton. 
Trong quá trình mô phỏng, phương pháp Monte – Carlo tập hợp tự hợp cần cập 
nhật phân bố của điện thế trong linh kiện ứng với một phân bố xác định của điện tích. 
Phân bố điện thế trong linh kiện có thể được xác định bằng việc giải phương trình 
Poisson, thông thường bằng phương pháp sai phân hữu hạn [4]. Khi đó việc giải phương 
trình Poisson chuyển thành việc giải một hệ phương trình tuyến tính cực lớn với hàng 
216
triệu phương trình và hàng triệu ẩn. Rõ ràng, việc giải hệ phương trình trên bằng một 
phương pháp giải tích là một việc bất khả thi và người ta phải sử dụng đến các phương 
pháp số. Đến nay nhiều phương pháp đã được xây dựng như: Jacobi, Gauss – Seidel, 
SOR, đa ô lưới (multigrid), iLU [5, 6]. Các phương pháp này có thể cho các kết quả 
chính xác tuy nhiên độ ổn định không cao và tốc độ hội tụ thấp. Gần đây, các phương 
pháp không gian con Krylov như CGS, BiCG, BICGSTAB, BICGSTAB2, 
BICGSTAB(l), GPBICG đã được phát triển và sử dụng như là các phương pháp hiệu 
quả trong việc giải các hệ phương trình tuyến tính thưa loại lớn [5]. Một số tác giả đã sử 
dụng phương pháp BICGSTAB để giải phương trình Poisson và đã thu được các kết quả 
chính xác với thời gian tính toán được rút ngắn nhiều lần [2, 7]. Đó là động lực để 
chúng tôi tiến hành tìm nghiệm của phương trình Poisson bằng phương pháp GPBICG, 
phương pháp được đánh giá là hoạt động ổn định hơn và cho kết quả nhanh hơn phương 
pháp BICGSTAB [8]. 
2. Giải phương trình Poisson ba chiều bằng thuật toán GPBICG 
Giả sử vật liệu là đồng nhất thì phương trình Poisson trong trường hợp ba chiều 
có dạng: 
2 2 2
2 2 2 ,
Sx y z

  
   
ở đây, là điện thế, là mật độ điện tích, S là hằng số điện môi tĩnh trong linh kiện; 
x , y , z là ba biến không gian. Để có thể dễ dàng thực hiện sai phân hữu hạn ta chia 
mô hình linh kiện thành các ô lưới và giả sử khoảng cách giữa các nút lưới theo các 
chiều không gian là bằng nhau, x y z . Tiến hành lấy sai phân hữu hạn phương 
trình (1) ta thu được hệ phương trình sau: 
, , 2
1, , , 1, , , 1 , , 1, , , 1, , , 16 ,
i j k
i j k i j k i j k i j k i j k i j k i j k
S
x
 
ở đây, 1, xi N , 1, yj N , 1, zk N với xN , yN , zN lần lượt là số nút lưới theo các 
chiều không gian Ox , Oy , Oz . Hệ phương trình (2) có thể được viết lại dưới dạng một 
phương trình ma trận như sau: ,A b (3) 
trong đó, ma trận A có dạng: 
   
     
     
   
1 1
2 2 2
Z 1 Z 1 Z 1
Z Z
,
0
0 N N N
N N
b c
a b c
A
a b c
a b

(1) 
(2) 
(4) 
217
với ja và jc là các ma trận một đường chéo chính: 
,
0
0
j ja c





còn jb là ma trận ba đường chéo chính: 
2(1 ) 1
1 2(1 ) 1
,
1 2(1 ) 1
1 2(1 )
0
0
jb





trong đó 2( ) 1x z . 
Bảng 1. Thuật toán GPBICG để tìm nghiệm của phương trình Poisson 
Đây là một phương trình ma trận loại lớn và việc giải phương trình này khá phức 
tạp. Dù dùng phương pháp nào thì để có thể giải hệ này với cách giải thông thường ta 
cũng đều cần một máy tính mạnh với bộ nhớ cực lớn để có thể lưu trữ và xử lý dữ liệu. 
(5) 
(6) 
218
May mắn là các phương pháp không gian con Krylov có thể hỗ trợ cách tính toán không 
cần lưu trữ các số liệu tính toán trung gian. Đây chính là ưu điểm lớn nhất của các 
phương pháp. Việc tính toán không cần dùng nhiều bộ nhớ có thể được thực hiện bằng 
cách khai triển các phép nhân ma trận thông qua hệ phương trình (4). Giải thuật 
GPBICG để tìm nghiệm của phương trình Poisson được khai triển trong Bảng 1. 
3. Kết quả mô phỏng và thảo luận 
Mô hình cấu trúc của đi-ốt p-i-n bán dẫn GaAs gồm một lớp bán dẫn thuần (i) 
kẹp giữa hai lớp bán dẫn pha tạp loại p và loại n như được chỉ ra trong Hình 1. Trong đó, 
mỗi lớp có độ dày tương ứng là id , pd và nd . Mật độ pha tạp acceptor và donor tương 
ứng là AN và DN , các tạp được phân bố từ bề mặt của các lớp p và n vào sâu bên trong 
linh kiện theo hàm phân bố Gauss. Trạng thái cân bằng nhiệt của linh kiện được xác lập 
bằng mô phỏng thời gian thực trước khi chiếu xung laser vào linh kiện. 
Hình 1. Mô hình đi-ốt p-i-n GaAs 
Chúng tôi đã sử dụng phương pháp Monte – Carlo tập hợp tự hợp để mô phỏng 
động lực học của hạt tải trong linh kiện trong trường hợp chiếu một xung laser với chiều 
dài của xung là 12 sf và năng lượng photon là 1.49 eV . Các tham số cấu trúc vùng 
năng lượng được sử dụng như sau: 1.42gapE eV
 , * 00.063em m , 
*
00.45hm m , 
*
00.222Lem m , và độ chêch lệch năng lượng giữa  và L 0.29LE eV . Chúng tôi 
giả thiết rằng 50p nd d nm , 340id nm và 
17 30.5 10AN cm
 , 
17 32.5 10DN cm
 và 16 35 10exN cm
 sau thời gian 1 ps . Kích thước theo ba chiều 
không gian của đi-ốt là 440 100 100x y zL L L nm nm nm , giả sử đi-ốt được nuôi 
cấy theo phương Ox . Mô hình linh kiện được chia thành các ô lưới không gian với 
khoảng cách giữa các nút lưới là 1050 10x y z m . Như vậy, ta sẽ có 89xN 
nút lưới theo phương Ox , 21yN nút lưới theo phương Oy và 21zN nút lưới theo 
phương Oz . Điện trường ngoài được đặt vào linh kiện dọc theo phương Ox và đi-ốt 
được phân cực nghịch, xem Hình 1. 
Hình 2 mô tả sự thay đổi vận tốc trôi dạt của điện tử theo các phương Ox , Oy 
và Oz và vận tốc trôi dạt toàn phần ứng với điện trường ngoài 100extE kV cm . Từ đồ 
219
thị ta thấy rằng, điện tử chủ yếu chuyển động trôi dạt theo phương Ox . Vận tốc trôi dạt 
toàn phần của điện tử được đóng góp chủ yếu từ thành phần vận tốc theo phương Ox , 
còn các thành phần vận tốc theo phương Oy và phương Oz cho đóng góp không đáng 
kể. Đặc biệt, tại thời điểm ban đầu sau khi chiếu xung laser vận tốc trôi dạt của điện tử 
theo phương Ox tăng nhanh vượt xa giá trị bão hòa rồi sau đó giảm nhanh về giá trị bão 
hòa, hiện tượng này được gọi là sự vượt quá vận tốc [1, 4]. 
Hình 2. Vận tốc trôi dạt của điện tử theo các phương khác nhau và vận tốc trôi dạt toàn phần 
như là hàm của thời gian ứng với 100extE kV cm 
Hình 3. Vận tốc trôi dạt của điện tử theo phương Ox như là hàm của thời gian ứng với các điện 
trường ngoài khác nhau 
220
Hình 3 mô tả sự phụ thuộc của vận tốc trôi dạt của điện tử theo thời gian ứng với 
các giá trị điện trường ngoài 70extE kV cm , 100kV cm và 130kV cm . Kết quả cho 
thấy, với điện trường ngoài càng cao thì sự vượt quá vận tốc xảy ra càng sớm và vận tốc 
càng nhanh chóng tiệm cận giá trị bão hòa. Trong cùng một khoảng thời gian, khi điện 
trường ngoài càng cao thì số điện tử nằm trong các trạng thái có thể tham gia vào quá 
trình tán xạ liên thung lũng càng lớn. Khi điện tử bị tán xạ từ thung lũng Γ sang thung 
lũng L vận tốc của điện tử bị giảm nhiều do khối lượng hiệu dụng của điện tử trong 
thung lũng L lớn hơn nhiều lần khối lượng hiệu dụng của điện tử trong thung lũng Γ. Hệ 
quả là vận tốc của điện tử càng giảm nhanh về giá trị bão hòa. 
Hình 4 cho kết quả so sánh vận tốc của điện tử thu được bằng hai chương trình 
mô phỏng ba chiều sử dụng thuật toán GPBICG và thuật toán BICGSTAB [2], cũng như 
kết quả thu được bằng chương trình mô phỏng một chiều [1]; 3 đồ thị gần như trùng 
nhau hoàn toàn. Điều đó cho thấy, phương pháp GPBICG hoạt động hiệu quả do 
chương trình một chiều đã được chứng minh là cho kết quả phù hợp với thực nghiệm. 
Hình 4. Vận tốc trôi dạt của điện tử theo phương Ox như là hàm của thời gian ứng với 
100extE kV cm : mô phỏng ba chiều (thuật toán GPBICG, BICGSTAB), và mô phỏng một 
chiều (thuật toán LU) 
Hình 5 mô tả sự phân bố điện thế không gian trong đi-ốt p-i-n bán dẫn GaAs 
theo hai phương Ox và Oy tại mặt cắt 10z nm ứng với điện trường ngoài 
100extE kV cm . Từ đồ thị ta thấy rằng điện thế trong linh kiện chủ yếu biến thiên theo 
phương Ox và gần như không đổi theo hai phương Oy và Oz . Kết quả này là hoàn 
toàn hợp lý do điện trường ngoài được đặt vào linh kiện theo phương Ox mà thôi. Kết 
quả này cũng phù hợp với các kết quả đã được công bố trước đây [1, 2]. 
221
Hình 5. Phân bố điện thế không gian trong đi-ốt p-i-n bán dẫn GaAs theo hai phương Ox và Oy 
tại mặt cắt 10z nm ứng với 100extE kV cm 
Hình 6. Sự phụ thuộc của chuẩn Euclid của vectơ thặng dư vào số vòng lặp của chương trình 
con Poisson ứng với 100extE kV cm 
Để so sánh tốc độ hội tụ và tính ổn định của thuật toán GPBICG so với thuật 
toán BICGSTAB chúng tôi đã tiến hành khảo sát sự phụ thuộc của chuẩn Euclid của 
vectơ thặng dư vào số vòng lặp của chương trình con Poisson, Hình 6. Chuẩn Euclid 
của vectơ thặng dư được tính theo công thức [5]. 
2
Tr r r , (7) 
với r b A là vectơ thặng dư. Từ đồ thị ta thấy rằng, thuật toán GPBICG cần ít vòng 
lặp hơn để tìm ra nghiệm có cùng chuẩn Euclid của vector thặng dư với thuật toán 
222
BICGSTAB, nghĩa là nó có tốc độ hội tụ nhanh hơn. Hơn thế nữa, đồ thị tương ứng với 
thuật toán GPBICG trơn hơn đồ thị tương ứng với thuật toán BICGSTAB, hàm ý rằng 
thuật toán GPBICG cho kết quả ổn định hơn thuật toán BICGSTAB. 
4. Kết luận 
Chúng tôi đã xây dựng thành công một chương trình giải phương trình Poisson 
ba chiều dựa trên thuật toán GPBICG dùng để tích hợp trong chương trình mô phỏng 
linh kiện na-nô bán dẫn bằng phương pháp Monte – Carlo tập hợp tự hợp. Để khảo sát 
các đặc trưng của phương pháp chúng tôi đã tiến hành mô phỏng động lực học ba chiều 
của hạt tải trong các đi-ốt p-i-n bán dẫn GaAs và so sánh với các kết quả mô phỏng đã 
được công bố trước đây. Các kết quả chỉ ra rằng, chương trình giải phương trình Poisson 
dựa trên thuật toán GPBICG không những có tốc độ hội tụ nhanh mà còn có tính ổn 
định cao hơn các chương trình từng được sử dụng. 
TÀI LIỆU THAM KHẢO 
[1]. D. N. Thao, S. Katayama, and K. Tomizawa, Numerical simulation of THz radiation by 
coherent LO phonons in GaAs p-i-n diodes under high electric fields, Journal of the 
Physical Society of Japan 73, (2004), 3177 – 3181. 
[2]. L. H. Linh, 3D simulation of carrier dynamics in GaAs p-i-n diodes by means of Monte 
- Carlo method, Master Thesis, Hue University’s College of Education, 2009. 
[3]. G. Klatt, B. Surrer, D. Stephan, O. Schubert, M. Fischer, J. Faist, A. Leitenstorfer, R. 
Huber, and T. Dekorsy, Photo-Dember terahertz emitter excited with an Er: fiber laser, 
Appl. Phys. Lett. 98, (2011), 021114 – 021114 - 3. 
[4]. K. Tomizawa, Numerical simulation of submicron semiconductor devices, Artech 
House, Boston London, 1993. 
[5]. H. A. Vorst, Iterative Krylov methods for large linear systems, Cambridge University, 
2003. 
[6]. A. Greenbaum, Iterative methods for solving linear systems, SIAM, Philadelphia, 1997. 
[7]. G. Speyer, D. Vasileska and S.M. Goodnick, Efficient Poisson equation solvers for 
large scale 3D simulations, Technical Proceedings of the 2001 International 
Conference on Modeling and Simulation of Microsystems, Nanotech 2001, Vol. 1, 
(2001), 23 - 26. 
[8]. Shao-Liang Zhang, GPBi-CG: Generalized product-type methods based on Bi-CG for 
solving nonsymmetric linear systems, SIAM J. Sci. Comput., Vol 18, No. 2, (1997), 537 
- 551. 
223
THREE-DIMENSIONAL SIMULATION OF NANO SEMICONDUCTOR 
DEVICES USING GPBICG ALGORITHM FOR THE SOLUTION 
OF THE POISSON'S EQUATION 
Dinh Nhu Thao, Duong Thi Diem My, 
Nguyen Chau Phuong Thi, Ngo Thanh Thuy 
 College of Pedagogy, Hue University 
SUMMARY 
The paper presents a way to build a three-dimensional Poisson solver based on 
GPBICG algorithm for integrating into a Monte Carlo simulation program of nano 
semiconductor devices. The program is used to simulate the three-dimensional carrier dynamics 
in GaAs p-i-n diodes. The results obtained are totally consistent with the published ones [1, 2]. 
The results show that GPBICG solver has not only faster convergence but also higher stability 
than the previous solvers [2]. 
Keywords: Simulation of semiconductor devices, 3D Poisson’s equation, GPBICG 
algorithm, Monte – Carlo method.