CC BY
27
Доклады БГУИР
2011 № 2 (56)
УДК 539.128.417
МОДЕЛИРОВАНИЕ РАССЕЯНИЯ ЭЛЕКТРОНОВ НА ИОНИЗИРОВАННОЙ ПРИМЕСИ В ПОЛУПРОВОДНИКАХ И ПОЛУПРОВОДНИКОВЫХ СТРУКТУРАХ МЕТОДОМ МОНТЕ-КАРЛО
Д.С. СПЕРАНСКИЙ, В.М. БОРЗДОВ, Д.В. ПОЗДНЯКОВ
Белорусский государственный университет пр. Независимости, 4а, 220030, Беларусь
Поступила в редакцию 9 декабря 2010
Рассмотрены основные подходы к описанию процесса рассеяния электронов на ионизированной примеси при численном моделировании электрофизических свойств полупроводников и полупроводниковых структур методом Монте-Карло. Показана возможность корректного описания процесса примесного рассеяния с помощью модели Ридли, в рамках которой предложен эффективный алгоритм нахождения полярного угла рассеяния 9. Для сравнения рассмотренных модельных приближений рассчитаны и проанализированы гистограммы угловых распределений и плотности вероятности для величины cos(6) для легированного объемного кремния при различных значениях концентрации примеси N¡ и энергии электрона Е.
Ключевые слова: метод Монте-Карло, рассеяние на ионизированной примеси, угол рассеяния 9, угловое распределение.
Введение
При адекватном численном моделировании электрофизических свойств полупроводников и полупроводниковых структур необходимо учитывать множество различных физических эффектов, влияющих на перенос в них носителей заряда. Одним из наиболее перспективных методов такого моделирования является метод Монте-Карло. Важнейшим преимуществом данного метода является возможность использования точных матричных элементов перехода носителя заряда при его взаимодействии с рассеивателем для любого механизма рассеяния в полупроводнике [1-3]. А это, как следствие, дает возможность строить достаточно сложные модели переноса, с помощью которых удается получать результаты, находящиеся в хорошем согласии с экспериментальными данными.
При моделировании кинетических процессов в полупроводниках методом Монте-Карло особого внимания заслуживает рассеяние на ионизированной примеси, поскольку оно существенно влияет на кинетические свойства сильнолегированных полупроводниковых областей. В частности, этот вид рассеяния будет играть весьма существенную роль при рассмотрении субмикронных полевых МОП-транзисторов из-за высокого уровня легирования истоковых и стоковых областей, а также подложки.
В настоящей работе рассмотрены основные подходы к описанию процесса рассеяния электронов на ионизированной примеси при численном моделировании методом Монте-Карло электрофизических свойств полупроводниковых структур, а также предложен новый эффективный алгоритм нахождения полярного угла рассеяния электрона 9 после его взаимодействия с заряженным центром для модели Ридли.
Модели примесного рассеяния
Рассматривая классическую задачу рассеяния электрона на кулоновском потенциале иона, полное поперечное сечение рассеяния можно представить в виде [3]:
с = лЬ2 = 2я|"а(г (9 9' de', (1)
где 9 - угол рассеяния, Ь - прицельный параметр, (9) - дифференциальное поперечное сечение рассеяния. Хорошо известно, что из-за расходимости величины ос! (9), при 9 ^ 0 интеграл в выражении (1) расходится, поэтому полное сечение рассеяния для кулоновского взаимодействия стремится к бесконечности, что является следствием бесконечно большого радиуса действия кулоновских сил [1].
На сегодняшний день при моделировании кинетических процессов в полупроводниках методом Монте-Карло чаще всего используются три основные подхода, позволяющие в той или иной мере снять расходимость интеграла в формуле (1): хорошо известные модели Кону-элл-Вайскопфа (CW) и Брукса-Херринга (ВН), а также менее распространенная модель исключения третьего тела или метод Ридли.
Приближение Конуэлл-Вайскопфа состоит в ограничении сверху допустимых значений прицельного параметра Ь половиной среднего расстояния между примесными центрами [1, 4, 5]:
Ьтах = 12 N )А (2)
ч-1
"I
где N - концентрация ионизированной примеси.
Подход Брукса-Херринга учитывает экранирование кулоновского потенциала двучас-тичного взаимодействия. В этой модели кулоновский потенциал заменяется потенциалом Юка-вы V(г) [1, 4, 5]:
Ze2
V (г) = ",-ехр(-фя), (3)
4Л880Г
( е2N ^
где Рх = -— - обратный радиус экранирования Дебая, г - расстояние до примесного
V 88о квТ )
центра, в0 - диэлектрическая постоянная, е - относительная диэлектрическая проницаемость полупроводника, Ze - заряд иона, N0 - концентрация свободных носителей, кв - постоянная Больцмана.
Однако следует иметь ввиду, что экранирование может быть не настолько сильным, чтобы ограничить величину поперечного сечения площадью круга с диаметром меньшим, чем среднее расстояние между примесями. Очевидно, что в этом случае процесс рассеяния уже некорректно рассматривать как простое двухчастичное взаимодействие [4]. Таким образом, хотя модель ВН строится, исходя из принципа двухчастичного рассеяния на ближайшем примесном центре, она все же не исключает возможности рассеяния электрона на более удаленном, третьем центре рассеяния. Если при этом ввести весовой множитель, который будет играть роль вероятности того, что данный акт рассеяния является чисто двухчастичным процессом рассеяния на ближайшем центре рассеяния, то подходы Брукса-Херринга и Конуэлл-Вайскопфа возможно согласовать [4]. Такой метод получил название эффекта исключения третьего тела или метода Ридли [1, 4, 6]. Его основная суть состоит в следующем.
Если предположить, что вероятность р существования какого-либо другого центра рассеяния, по отношению к которому прицельный параметр Ь будет находиться в интервале [Ь,Ь+ёЬ], равна [4]:
р = 2жNIaЬdЬ , (4)
то несложно показать [4], что вероятность отсутствия рассеивающих центров с прицельным параметром, меньшим Ьт, можно представить в виде:
Р(Ь) = ехр(-каЬ2т N),
(5)
где а - среднее расстояние между ионами.
Тогда, если дифференциальное поперечное сечение, рассчитанное по модели ВН, умножить на эту вероятность Р(Ь), то приходим к простой формуле для определения соответствующего дифференциального сечения по модели Ридли [6]:
^ (к, е)=Р (ь (к, е).
(6)
Таким образом, вероятность Р(Ь) является множителем, исключающим влияние третьего тела.
На рис. 1 приведены зависимости интенсивностей примесного рассеяния в кремнии от энергии электрона Е при различных концентрациях N и температурах Т. Если учесть, что ку-лоновское рассеяние является остроугловым и требует больших затрат машинного времени [1], то в этой связи использование в алгоритме Монте-Карло модели ВН становится более предпочтительным по сравнению с моделью CW в широком диапазоне энергий электрона, т.к. в первом случае вероятность рассеяния уменьшается при достаточно высоких энергиях. Однако с уменьшением концентрации примеси модель ВН дает неоправданный рост интенсивности рассеяния при малых значениях энергии. Это, очевидно, связано с тем, что использование первого борновского приближения, в условиях которого она получена, становится некорректным. С учетом вышесказанного более приемлемым в данном случае является использование модели Ридли, которая приводит к результатам, находящимся в лучшем согласии с экспериментальными данными [6].
10»
№
ю1-
ш11
■ -г* 0 г
■ -<Г тЦоок
- Т=77К
:-,т1= 1.0
< >
--1-1—1 11111 -1— - — —►
Е.эВ
(С № ю-3 10"! 1
Рис. 1. Зависимость интенсивностей рассеяния на ионизированной примеси от энергии носителя заряда в кремнии при двух значениях N и Т; Z=1: (—)модель Конуэлл Вайскопфа (CW), (-) модель Брукса-Херринга (ВН), (.......) модель Ридли [6]
Моделирование угловых распределений
Важную роль при моделировании методом Монте-Карло взаимодействия электрона с ионом примеси играет выбор состояния электрона после рассеяния. Рассматривая этот процесс абсолютно упругим, оно задается единственной величиной - полярным углом 0, а азимутальный угол ф считается равномерно распределенной случайной величиной на интервале от 0 до 2п [1, 5, 8]. Ниже кратко рассмотрим процедуры розыгрыша угла 0.
Как известно, интенсивность рассеяния на примеси в первом борновском приближении для потенциала Юкавы (или вероятность рассеяния в единицу времени) можно представить в виде [8]:
Ре,г (Е) =
22ЫтеАт * k г^е,
2лЯг в„
^соэ -1
d соэе'
[р2 + 2k2(1 - соэ е')]2'
(7)
где т* - эффективная масса электрона, Н - редуцированная постоянная Планка.
Используя метод обратных функций, несложно получить следующее выражение для розыгрыша методом Монте-Карло угла рассеяния 0 [8]:
I
соэ е, соэ е.
dсоэе'
г =
Р (Е)_ е [р 2 + 2k2 (1 - соэ е')
Ре1 (Е) I - е
d соэе'
(8)
[р 2 + 2k2 (1 - соэ е')]
где г - равномерно распределенное на интервале [0, 1] случайное число. Из (8) непосредственно следует:
соэ 0 = 1 - 2
(1 + 4 ЕЕ ) эт^е^А) + г соэ2 (етп/2)
1 - 2А_/ Ер>_
1 + 4(%в )[1 - гсоэ2 (етп/2)]
(9)
где Еп =
2т *
Для модели ВН, учитывающей экранирование, полагая 0ш1п = 0, выражение (9) примет вид [1, 8]:
соэ е^я = 1 - [2(1 - г)/(4 (Е/Еп ) г +1)] .
Для модели CW Рх = 0 и тогда [8]:
соэ е™ =
где Еь =
((Е/Еь )2 г -1)/((Е/Еь )2 г +1)
(10)
(11)
8^0Ьтах
Для модели Ридли, с учетом выражения (5), несложно получить следующее выражение для расчета вероятности рассеяния в единицу времени [1, 7]:
Ря (Е) = -а
2Е
т'
(
1 - ехр
ау]т
<(Е)
42Е
(12)
где Рвн (Е) - вероятность рассеяния в единицу времени, рассчитанная по модели Брукса-Херринга.
В этом случае в работе [7] для розыгрыша угла рассеяния 0 была предложена процедура, согласно которой сначала разыгрывается прицельный параметр Ь, а затем путем инверсии выражения (1), где ака (е) = овНс1 определяется значение 0. Согласно [7], ис-
пользуя метод обратных функций, будем иметь:
I ь Ь ехр^ъаЫ^^Ь
■ * ь /
?ьт 2 ■
I Ь ехр(-%аЫ1Ь ^Ь
2
2
е
Поскольку
bm 2 1 2
I b exp(-%aNIb )db =-(1 - exp(-naN¡bm)),
J° 2naN I m
то [7]:
z = sxp{-%aNjb2) - exp(-TCaN7b,m) (14)
1 - exp(-%aN¡b,)
где bm определяется из соотношения:
2 f1
nbm = 2nJ_i oK,d (cos(0))d cos(0).
Из (14) следует, что параметр b, очевидно, может быть найден с помощью метода Монте-Карло.
Для оптимизации процедуры нахождения угла 0 в методе Ридли, чтобы избежать вначале розыгрыша параметра b согласно (14), а затем процедуры инверсии уравнения (1), в настоящей статье предлагается непосредственно определять угол из выражения (9) с учетом поправки (5), учитывающей вероятность отсутствия другого центра рассеяния. Тогда, с использованием выражения (5) , на основании (7) и (13) имеем:
•cos 0™п d cos 0'
|c 2
cos0z [P2 + 2k2 (1 - cos0')] ^ exp(-%aN¡b2) - exp(-naN^
, (15)
1 2 rcos 0min_d cos 0'__1 - exp(-naNIb2m) '
J-1 [P2 + 2k2 (1 - cos 0')]2
где z1 и z2 два случайных числа, равномерно распределенных на отрезке от 0 до 1. Таким образом, формула (15) позволяет непосредственно определить угол 0, разыгрывая всякий раз два случайных числа z1 и z2. Эта процедура значительно упрощается, если в правой части (15) экспоненты разложить в ряд Тейлора и ограничиться двумя членами этого ряда. В этом случае можно получить аналитическое выражение для cos (0). Действительно,
для модели Ридли можно использовать простую связь прицельного параметра b и угла рассеяния 0 [1]:
b = W1 + cos(0) , (16)
л/1 - cos(0)
Ze2
где K =-- . Выполняя операции интегрирования в выражении (15), а также учитывая
4nss°m * v
(16), приходим к следующему уравнению: (cos0г -1)2(1 + cos0m)(2EzZ2(1 + cos0m) - (E - ED)(1 - cos0m)) -(cos 0г -1) [Z1 z2Ed (1 + cos 0m)2 + 4(Ed + 4E)(1 - cos 0m)] + 2(Ed + 4E)(1 - cos2 0m) = 0.
Введем обозначения: m = (1 + cos 0m )(2Ez z2(1 + cos 0m) - (E - Ed )(1 - cos 0m), n = Z1 z2Ed (1 + cos 0m )2 + 4(Ed + 4E)(1 - cos 0m), h = 2(Ed + 4E)(1 - cos2 0m).
(17)
n — \J n2 — 4mh
Решением уравнения (17) очевидно будет корень вида -. Второй корень
2m
n + \J n2 — 4mh „
- не имеет смысла, поскольку в этом случае прицельный параметр b больше мак-
2m
симально допустимого значения.
Таким образом, окончательно получаем аналитическое выражение для розыгрыша cos(0z):
л Vn2 — 4mh — n
cos 0Z = 1 н--. (18)
2m
Результаты расчетов угла рассеяния 8
Для подтверждения корректности предложенной процедуры розыгрыша угла 0, в настоящей работе в качестве примера были рассчитаны плотности вероятности ^Р(соб (6)) при
Т = 300К, ^=1024м_3 и двух различных значениях энергии Е (рис. 2, 3).
Р(СО5|0)>
cos|0J
Общая модель [8]- - -Наш подход Модель Ридли согласно[7]
Рис. 2. Плотность вероятности cos(6) при £=0,1 эВ, NI= 1024 м-3, 7=300 K
P(COS|0)|
/
/1
<У/
—............... 1 |П ....................|1
СИ СО К CD 1Л ^ П h- СО СГ|
срсрсрсрсрсрсрсрср ооооооооо
COS(0)
Общая модель [8]- - -Наш подход —Модель Ридли согласно[7]
Рис. 3. Плотность вероятности соз(9) при Е=1 эВ, Ж=1024 м-3, Т=300 К
Из этих рисунков видно, что предложенная модель хорошо согласуется с результатами, рассчитанными по общей формуле [8] и по процедуре розыгрыша угла 0, предложенной в [7].
Угловое распределение, полученное нами согласно формуле (15) при энергии электрона Е=0,01 эВ (рис. 4) носит выраженный остроугловой характер, что говорит об адекватном описании процесса рассеяния электронов на ионизированной примеси [1].
70 60 50 40 30 20 10 0
-1 -0 9 -0,8 -0.7 -0,6 -0,5 -0,4 -0,3 -0,2 -0 1 0 0 1 0.2 0.3 0.4 0,5 0,6 0,7 0,8 0,9 1
COS(0)
■ Модель ВН *Уодегь CW ■ Наш подход
Рис. 4. Гистограмма распределения cos6r при N=1023 м-3, £=0,01 эВ, T=300 K
Отметим, что предложенная в настоящей статье процедура розыгрыша углового распределения является более общей, чем та, которая была описана в [9].
Заключение
Таким образом, предложенный в данной статье подход позволяет упростить процедуру розыгрыша углового распределения при моделировании процессов рассеяния электронов на ионизированной примеси методом Монте-Карло по модели Ридли. Такой подход позволяет непосредственно определять угол рассеяния из общей формулы (9), минуя процедуру инверсии, описанную в работах [7] и [10].
MONTE-CARLO SIMULATION OF IONIZED IMPURITY SCATTERING IN SEMICONDUCTORS AND SEMICONDUCTOR STRUCTURES
D.S. SPERANSKY, V.M. BORZDOV, D.V. POZDNYAKOV
Abstract
The basic approaches to the description of ionized impurity scattering process at numerical simulation of electrophysical properties of semiconductors and semiconductor structures by Monte-Carlo method are considered. Possibility of its correct description within the frame of Ridley model is shown. Effective algorithm of polar angle 9 selection in case of Ridley model is offered.
As an example doped bulk silicon angle distribution histograms and probability density of cos(0) for various values of impurity concentration Nj and electrons energy E are presented.
Литература
1. Иващенко В.М., Митин В.В. Моделирование кинетических явлений в полупроводниках. Метод Монте-Карло. Киев, 1990.
2. Борздов В.М., Комаров Ф.Ф. Моделирование электрофизических свойств твердотельных слоистых структур интегральной электроники. Минск, 1999.
3. Борздов В.М. Моделирование методом Монте-Карло приборных структур интегральной электроники. Минск, 2007.
4. Ридли Б. Квантовые процессы в полупроводниках. М., 1986.
5. Jacoboni C., Lugly P. The Monte Carlo Method for Semiconductor Device Simulation. Wien-New York, 1989.
6. Speransky D.S. Monte Carlo simulation of ionized impurity scattering process in bulk silicon. NDTCS'08
7. Van deRoer T.G., Widdershoven F.R. // J. Appl. Phys. 1986. V. 59, №3. P. 813-815.
8. Jacoboni C., Reggiani L. // Rev. Mod. Phys. 1983. V. 55, №3. P. 645-705.
9. Борздов В.М., Сперанский Д.С. // CriMiCo'2010. 2010. P. 867-868.
10. J. Bude // Monte Carlo Device Simulation: Full Band and Beyong, edited by Karl Hess, Kluwer Academic Publishers, Boston/Dordrecht/London
Г Г Г1 ' 'I ГГГ V I 'I 'I' I l"l Y V Г1 V Y Г1 Yl"l
"i11! V "V'm'M'M?M
ПГПГ/
