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

TẠP CHÍ KHOA HỌC, Đại học Huế, Số 65, 2011<br /> MÔ PHỎNG BA CHIỀU LINH KIỆN NA-NÔ BÁN DẪN<br /> VỚI LỜI GIẢI PHƯƠNG TRÌNH POISSON<br /> DỰA TRÊN THUẬT TOÁN GPBICG<br /> Đinh Như Thảo, Dương Thị Diễm My,<br /> Nguyễn Châu Phương Thi, Ngô Thanh Thủy<br /> Trường Đại học Sư phạm, Đại học Huế<br /> <br /> TÓM TẮT<br /> 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<br /> 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<br /> 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ô<br /> 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ô<br /> 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<br /> đâ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<br /> 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<br /> từng được sử dụng [2].<br /> 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<br /> GPBICG, phương pháp Monte – Carlo.<br /> <br /> 1. Giới thiệu<br /> 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<br /> 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<br /> 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ệ<br /> 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<br /> 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<br /> 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<br /> 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<br /> 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<br /> 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<br /> 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.<br /> Trong quá trình mô phỏng, phương pháp Monte – Carlo tập hợp tự hợp cần cập<br /> 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.<br /> 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<br /> 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<br /> 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<br /> 215<br /> <br /> 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<br /> 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<br /> pháp số. Đến nay nhiều phương pháp đã được xây dựng như: Jacobi, Gauss – Seidel,<br /> SOR, đa ô lưới (multigrid), iLU [5, 6]. Các phương pháp này có thể cho các kết quả<br /> 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<br /> pháp không gian con Krylov như CGS, BiCG, BICGSTAB, BICGSTAB2,<br /> BICGSTAB(l), GPBICG đã được phát triển và sử dụng như là các phương pháp hiệu<br /> 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ử<br /> dụng phương pháp BICGSTAB để giải phương trình Poisson và đã thu được các kết quả<br /> 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 để<br /> 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,<br /> 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<br /> pháp BICGSTAB [8].<br /> 2. Giải phương trình Poisson ba chiều bằng thuật toán GPBICG<br /> Giả sử vật liệu là đồng nhất thì phương trình Poisson trong trường hợp ba chiều<br /> có dạng:<br /> <br />  2  2  2<br /> <br />  2  2  ,<br /> 2<br /> x<br /> y<br /> z<br /> S<br /> <br /> (1)<br /> <br /> ở đâ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;<br /> 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<br /> 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<br /> 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<br /> trình (1) ta thu được hệ phương trình sau:<br /> <br /> i 1, j ,k  i , j 1,k  i , j ,k 1  6i , j ,k  i 1, j ,k  i , j 1,k  i , j ,k 1  <br /> <br /> i , j ,k 2<br /> x , (2)<br /> S<br /> <br /> ở đây, i  1, N x , j  1, N y , k  1, N z với N x , N y , N z lần lượt là số nút lưới theo các<br /> 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<br /> phương trình ma trận như sau:<br /> A  b, (3)<br /> trong đó, ma trận A có dạng:<br />  b1 <br /> <br />  a2 <br /> A<br /> <br /> <br /> <br /> <br /> <br />  c1 <br /> b2   c2 <br /> <br /> 0<br /> <br /> <br /> 0<br /> <br />  aNZ1  bNZ1 <br />  a NZ <br /> 216<br /> <br /> <br /> <br /> <br /> ,<br /> <br /> cNZ1 <br /> bNZ  <br /> <br /> (4)<br /> <br /> với  a j  và  c j  là các ma trận một đường chéo chính:<br /> <br /> <br /> <br /> <br />  a j   c j   <br /> <br /> <br /> <br /> <br /> <br /> 0<br /> <br /> <br /> <br /> <br /> 0<br /> <br /> <br /> <br /> <br /> <br /> <br /> ,<br /> <br /> <br /> <br />  <br /> <br /> (5)<br /> <br /> còn b j  là ma trận ba đường chéo chính:<br /> <br /> 1<br />  2(1   )<br /> <br /> <br /> <br /> 1<br /> 2(1   ) 1<br /> <br /> <br /> ,<br /> b j   <br /> <br /> <br /> <br /> <br /> <br /> 1 2(1   )<br /> 1<br /> <br /> <br /> 1<br /> 2(1   ) <br /> <br /> <br /> 0<br /> <br /> (6)<br /> <br /> 0<br /> <br /> trong đó   (x z )2  1 .<br /> Bảng 1. Thuật toán GPBICG để tìm nghiệm của phương trình Poisson<br /> <br /> Đâ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<br /> 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<br /> 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.<br /> 217<br /> <br /> 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<br /> 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<br /> 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<br /> 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<br /> GPBICG để tìm nghiệm của phương trình Poisson được khai triển trong Bảng 1.<br /> 3. Kết quả mô phỏng và thảo luận<br /> 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)<br /> 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 đó,<br /> mỗi lớp có độ dày tương ứng là di , d p và d n . Mật độ pha tạp acceptor và donor tương<br /> ứng là N A và N D , 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<br /> 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<br /> bằng mô phỏng thời gian thực trước khi chiếu xung laser vào linh kiện.<br /> <br /> Hình 1. Mô hình đi-ốt p-i-n GaAs<br /> <br /> Chúng tôi đã sử dụng phương pháp Monte – Carlo tập hợp tự hợp để mô phỏng<br /> độ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<br /> dài của xung là 12 fs và năng lượng photon là 1.49 eV . Các tham số cấu trúc vùng<br /> <br /> năng lượng được sử dụng như sau: E gap<br />  1.42 eV , m* e  0.063 m0 , m* h  0.45 m0 ,<br /> <br /> *<br /> mLe<br />  0.222 m0 , và độ chêch lệch năng lượng giữa  và L EL  0.29 eV . Chúng tôi<br /> <br /> giả<br /> <br /> thiết<br /> <br /> rằng<br /> <br /> d p  d n  50 nm<br /> <br /> ,<br /> <br /> di  340 nm<br /> <br /> và<br /> <br /> N A  0.5  1017 cm 3<br /> <br /> ,<br /> <br /> N D  2.5  1017 cm3 và N ex  5  1016 cm 3 sau thời gian 1 ps . Kích thước theo ba chiều<br /> <br /> không gian của đi-ốt là Lx  Ly  Lz  440 nm  100 nm  100 nm , giả sử đi-ốt được nuôi<br /> 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<br /> khoảng cách giữa các nút lưới là x  y  z  50 10 10 m . Như vậy, ta sẽ có N x  89<br /> nút lưới theo phương Ox , N y  21 nút lưới theo phương Oy và N z  21 nút lưới theo<br /> 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<br /> được phân cực nghịch, xem Hình 1.<br /> 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<br /> và Oz và vận tốc trôi dạt toàn phần ứng với điện trường ngoài Eext  100 kV cm . Từ đồ<br /> 218<br /> <br /> 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<br /> 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 ,<br /> 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<br /> 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ử<br /> 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<br /> hòa, hiện tượng này được gọi là sự vượt quá vận tốc [1, 4].<br /> <br /> 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<br /> như là hàm của thời gian ứng với Eext  100 kV cm<br /> <br /> 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<br /> trường ngoài khác nhau<br /> 219<br /> <br />