CC BY
49
Вычислительные технологии
Том 18, № 5, 2013
Численная модель ВЧ-разряда
>к
в плазмохимическом реакторе*
Ю. Н. Григорьев, А. Г. Горовчук Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: grigor@ict.nsc.ru, alg@eml.ru
Рассматривается численная модель аксиально-симметричного высокочастотного разряда в гидродинамическом приближении. Для решения уравнений непрерывности для электронов и ионов и уравнения энергии электронов использована неявная экспоненциальная разностная схема, обеспечивающая положительные значения электронной температуры и концентраций плазменных компонентов. Разработанный алгоритм расчёта ВЧ-разряда замыкает численную модель плазмохими-ческого травления, созданную в предыдущих работах авторов.
Ключевые слова: высокочастотный разряд, гидродинамическая модель, плаз-мохимическое травление.
Введение
В производстве изделий микроэлектроники важнейшим инструментом являются плаз-мохимические реакторы (ПХР), процесс травления в которых осуществляется в плазме высокочастотного (ВЧ) разряда. Электрические характеристики плазмы, в частности концентрации и энергия электронов и ионов, определяются такими факторами, как конструкция реактора, вкладываемая в разряд мощность, его частота и др. Изменение концентраций электронов и ионов в разряде может существенно повлиять на скорость и однородность травления образцов. В связи с новыми требованиями к качеству обработки микросхем в настоящее время широко развивается исследование высокочастотного ВЧ-разряда применительно к условиям ПХР. Особую значимость приобретают численные расчёты ВЧ-разряда в двумерной постановке, учитывающие изменение электронной плотности вдоль электрода, несущего обрабатываемые образцы.
Исследование физических процессов в ВЧ-разряде — сложная многопараметрическая задача. В оптимизационных расчётах реальных конструкций ПХР обычно используются экспериментальные данные по концентрациям заряженных частиц. Простые модели процессов травления основаны на аналитических приближениях, позволяющих на качественном уровне описать основные особенности структуры ВЧ-разряда [1]. Более реалистичные подходы к расчёту разряда, как правило, представляют собой самостоятельные задачи и до настоящего времени носят преимущественно исследовательский характер.
Детальное описание плазменных процессов в ВЧ-разряде даёт решение кинетических уравнений Больцмана для функций распределения нейтральных частиц, ионов и электронов, что позволяет определить характеристики плазмы в широком диапазоне
* Работа выполнена при финансовой поддержке РФФИ (проект № 11-01-00064) и гранта Президента РФ для государственной поддержки ведущей научной школы РФ № НШ-6293.2012.9.
температур и давлений. Однако отсутствие достоверных данных по сечениям основных кинетических процессов и высокие требования к вычислительным ресурсам делают данный подход не пригодным для практических целей. Эти трудности частично преодолеваются путём использования гибридных моделей, объединяющих гидродинамический и кинетический подходы [2], в которых функция распределения электронов по энергиям находится из уравнения Больцмана, а концентрации ионов, электронов и нейтральных частиц вычисляются из соответствующих уравнений непрерывности.
Относительно полное представление о распределении концентраций электронов и ионов, а также температуры электронов даёт моделирование ВЧ-разряда в гидродинамическом приближении [3, 4]. При этом рассматриваются уравнения переноса для гидродинамических моментов — плотности, скорости и энергии заряженных частиц, полученные интегрированием кинетических уравнений Больцмана для ионов и электронов, и уравнение Пуассона для самосогласованного электрического потенциала. Гидродинамические модели применяются в диапазоне рабочих давлений 0.1 ^ 3.0 Торр, при котором средняя длина свободного пробега нейтральных частиц, ионов и электронов значительно меньше толщины приэлектродных слоёв.
К преимуществам гидродинамического подхода можно отнести использование конечно-разностных схем, требующих сравнительно небольшого времени расчёта до выхода ВЧ-разряда на периодический режим, особенно при учёте большого количества процессов взаимодействия частиц. Однако временной масштаб некоторых реакций может быть достаточно велик, что требует интегрирования основных уравнений до ~ 105 циклов ВЧ-разряда [4].Тем не менее умеренная вычислительная трудоёмкость гидродинамического подхода позволяет применгть его для адекватного исследования реакторных процессов.
Анализ устойчивости конечно-разностных схем, используемых в расчётах по гидродинамической модели, даёт различные ограничения на шаг по времени [5], что определяет области их применимости. На стадии зажигания разряда при малой проводимости плазмы и небольших градиентах электрического поля для решения уравнений непрерывности можно использовать явные конечно-разностные схемы первого порядка точности по времени и пространству, обеспечивающие удовлетворительные результаты при числе Куранта для электронов Ке = veт/h < 1, где ve — дрейфовая скорость электронов, т, h — шаги по времени и пространству [5]. Основные уравнения обычно решаются последовательно без итераций самосогласования.
На стадии горения, характеризующейся образованием тонких приэлектродных сло-ёв и квазинейтрального столба, из-за сильных градиентов параметров плазмы и большой проводимости проведение расчётов по явным схемам становится крайне неэкономичным. Вследствие ограничений, накладываемых условием Куранта, шаг по времени оказывается много меньше характерного времени развития разряда, что может потребовать более 104 — 105 шагов по времени до выхода на периодический режим. Здесь эффективно применение неявных конечно-разностных схем, позволяющих использовать шаг по времени в ~ 10 — 102 раз больший, чем даёт условие Куранта. Однако в численных схемах, в которых распределение потенциала находится непосредственно из уравнения Пуассона, возникает дополнительное ограничение на максимальный шаг по времени вида т < 1/4по, где о — проводимость плазмы [5, 6]. Данное ограничение даёт величину временного шага на порядок больше, чем условие Куранта. Кроме того, само решение системы нелинейных разностных уравнений на каждом шаге по времени в этом случае сопряжено со значительными трудностями.
В настоящее время при моделировании газовой плазмы напряжённость электрического поля вычисляется из уравнения сохранения электрического тока. Для решения данного уравнения применяются консервативные по заряду разностные схемы, в которых шаг по времени ограничивается периодом лэнгмюровских колебаний ионов [6]. При использовании чисто неявных схем для уравнения тока решение разностных уравнений требует дополнительных итераций по самосогласованию, что снижает эффективность вычислительного алгоритма. Применение в нём явных численных схем приводит к известному ограничению на шаг по времени. Поэтому для построения эффективных численных алгоритмов рациональным подходом является использование преимуществ каждой из схем. Так, в [5] предложена безытерационная полунеявная схема, где в области квазинейтрального столба плотности электронов и ионов вычисляются по явной схеме, а в оставшихся точках — по неявной. Сравнение показывает, что расчёты по полунеявной схеме в несколько раз быстрее, чем по чисто неявной. Однако в узлах сетки при переходе от явного вычисления концентраций частиц к неявному схема обладает условной аппроксимацией.
Проведённый анализ позволил определить структуру эффективного алгоритма численного моделирования ВЧ-разряда в гидродинамическом приближении. В настоящей работе рассматривается безытерационный численный алгоритм, в котором для решения уравнений непрерывности электронов и ионов и уравнения баланса электронной энергии используется неявная экспоненциальная конечно-разностная схема. Уравнение Пуассона решается с применением конечно-разностной схемы типа "крест", к сожалению, несколько ограничивающей эффективность численного алгоритма. Для ускорения расчётов выполнено распараллеливание вычислительного алгоритма на кластерной системе с применением MPI. В расчётах на основе предлагаемого подхода исследуются характеристики плазмы ВЧ-разряда для плоской и цилиндрической конфигурации электродов и их влияние на реакторные процессы травления.
1. Математическая модель
Численно моделируется аксиально-симметричный высокочастотный разряд в реакторе плазмохимического травления. В расчётах используется гидродинамическая модель аксиально-симметричного ВЧ-разряда, включающая уравнения непрерывности для электронов и положительных ионов, уравнение баланса энергии электронов и уравнение Пуассона для электрического потенциала [3, 4]. Начально-краевая задача рассматривается в цилиндрической области, где электродами являются торцевые поверхности реакционной камеры.
Вводятся следующие переменные:
Щ VI ~ А
Рг = —, VI = —, А = —, I = е,р,
пе0 V10 А10
, Те ф г г Ь
Те 0 фо П Ь п
Здесь индексы е,р означают электроны и ионы, щ — плотность частиц; VI, А — подвижность и коэффициент диффузии; Те — электронная температура; ф — потенциал; / — частота активации; £ — время. Плотность, температура, подвижность, коэффициент диффузии и потенциал обезразмеривались на среднюю плотность электронов пе0,
температуру Te0, соответствующую тепловой энергии электронов, характерные значения , Di0 и потенциал ВЧ электрода ф0. Цилиндрические координаты r, z нормированы на радиус цилиндра rt и межэлектродное расстояние L.
Концентрации плазменных компонентов находятся из уравнений непрерывности для электронов и ионов
.dpi , Л2(— , -М ,
где
+A4#+V+*=^ (1)
T = ñ = í 1, l = e,
De0 , Pl 1 De/Dp, l = p. Потоки частиц определяются формулами
dDipi dDipi ¡ -1, l = e, . ,
ji€ = -siPeiAiPi- "^F", Лс = -siPei^iPidz - , s = \ 1, l = p. (2)
Нижние индексы £ и Z в обозначениях потоков относятся к направлениям вдоль радиальной и аксиальной координат. Число Пекле характеризует соотношение процессов переноса и диффузии:
Pei = ~ñ~. Dio
Источниковый член описывает ионизацию исходных молекул и процесс прилипания электронов к ионам:
Pe (Dar i - Га), l = e, u . . . 1 va¿2
S = | peDar.; / = p',r' = exp[(1 - ^ E'/kT'0], Deo •
где r. — скорость ионизации, ra — скорость прилипания, E. — энергия ионизации, k — постоянная Больцмана, va — частота прилипания. Число Дамкелера, пропорциональное отношению интенсивности ионизации к диффузии, вычисляется по формуле
L2
Da = —-nk.o exp (-E./kTeo),
De0
здесь n — объёмная плотность газа, k;o — коэффициент скорости ионизации. Распределение электрического потенциала находится из уравнения Пуассона
21 5 (i
Введённый параметр Г характеризует напряжённость электрического поля:
L2 eñeü
г
где е — заряд электрона, е — диэлектрическая постоянная в вакууме.
Электронная температура вычисляется из уравнения баланса электронной энергии
3 Т ддТ <*'•) + А2 (^ + |) + ^ - X ( | + |) + ад = 0, (4)
где компоненты потока электронной энергии определяются соотношениями
3% = 3*е€ йе, 3-9С = Зес «е.
Параметр х характеризует отношение энергии электрического поля к тепловой энергии электронов, параметр в отвечает за потери энергии в реакциях ударной ионизации. Их выражения имеют вид
х
ер0 к
е0
в = ^, ке0
к
е0
2 кТе 0.
Здесь к — потери энергии при ионизации, ке0 — тепловая энергия электронов.
Низкотемпературная плазма находится между высокочастотными электродами, и электронная и ионная плотности быстро снижаются вне радиуса этой зоны. Такая постановка соответствует ВЧ-разряду между двумя плоскими поверхностями. В действительности ВЧ-разряд инициируется в цилиндрическом объёме, в котором низкотемпературная плазма ограничивается непроницаемой боковой стенкой в радиальном направлении. В настоящей работе рассматриваются две структуры ВЧ-разряда — между двух плоских поверхностей и в цилиндрическом объёме. В первом случае разряд рассчитывается в области, ограниченной электродами (0 < С < С0, 0 < С < 1), где Се = г0/п — радиус электрода, во втором — во всем объёме реакционной камеры (0 < С < 1, 0 < £ < 1).
На границах расчётной области для уравнений (1)-(4) ставятся следующие краевые условия:
— на заземленном аноде (0 < С < С0, С = 0):
Ре = 0, дАРр = 0, р = 0, « = «е.;
дС
— на силовом катоде (0 < С < С0, С = 1):
дБ р рр
Зе^ ,
д(
0, р = вШ (2ПТ) , $е = дез]
— на оси симметрии (С = 0, 0 < £ < 1)
дде
дРе =0 дРр = 0 дР = 0 дде =0-дС 0, дС 0 дС 0, дС 0;
для ВЧ-разряда между плоскостями:
— на выходе из плазменной области (С = С0, 0 <С< 1):
дРе =0 дБ рРр =0 д^ = 0 дде = 0
дС 0, дС 0, дС 0 дС 0;
для ВЧ-разряда в цилиндрическом объёме:
— на боковой поверхности (С =1, 0 < С < 1):
Ре = 0, ^ = 0, | = 0, « = «е.!
в верхней части реакционной камеры (£0 < £ < 1, £ =1):
Ре = 0, ^ = 0, ^ = 0, = ^ на перфорированной границе откачки (£0 < £ < 1, С = 0):
Ре = 0, ^ = 0, ^ = 0, ^ = ^
Здесь 7 — коэффициент вторичной электронной эмиссии, — электронная температура на электродах.
В силу быстрой рекомбинации электронов на проводящей поверхности анода электронная плотность на нём полагалась равной нулю. На катоде учитывался слабый поток вторичной электронной эмиссии. Краевые условия для ионов записывались в виде диффузионных потоков к границам плазменной области. На оси реактора задавались условия симметрии. Потенциалы и электронная температура на электродах предполагались заданными. При моделировании ВЧ-разряда между двумя плоскими поверхностями на границе разрядной области по радиусу принимались мягкие краевые условия для электронной и ионной плотности и электронной температуры. При рассмотрении ВЧ-разряда в цилиндрическом объёме на боковой стенке реактора в верхней части реакционной камеры и на перфорированной границе откачки плотность электронов полагалась равной нулю, электронная температура — равной температуре на электродах, а для потенциала использовались мягкие краевые условия.
Численное решение уравнений (1)-(4) является трудной вычислительной задачей из-за сильной нелинейности и взаимосвязи уравнений. Уравнения непрерывности (1) содержат большие числа Пекле (до 104), что требует применения специальных численных методов. Интегрирование уравнений с малым параметром при старшей производной осложняется образованием пограничных слоёв с большими градиентами искомых величин, ведущими к неравномерной сходимости решения. В качестве одного из подходов к решению возникающей проблемы используются повышение порядка точности разностной схемы и подробные сетки в области сильного изменения решения. Однако с увеличением жёсткости уравнения эффективность такого подхода существенно снижается. Кроме того, раздельная аппроксимация конвективных и диффузионных членов центральными или направленными против потока разностями для первых производных и центральными разностями для вторых производных практически непригодна при построении монотонных схем.
Основным методологическим подходом при построении безусловно монотонных разностных схем является их регуляризация, частный случай которой — метод экспоненциальной подгонки или интегральных тождеств [7-9] — используется в настоящей работе. Для численного решения уравнений (1), (4), записанных для суммарных потоков, включающих конвекцию и диффузию, применяется неявная экспоненциальная разностная схема, предложенная для одномерного случая в [10, 11]. На рис. 1 показан шаблон использованной сетки. Чёрные квадраты соответствуют узлам с целыми индексами, белые и чёрные кружки — узлам с полуцелыми индексами на радиальной и осевой координатах. Грани ячеек (шаги) сетки с целыми индексами обозначены как Д^+1/2 = £¿+1 — £ и Д^.+1/2 = 0+1 — . Для каждого типа узлов приведены обозначения вычисляемых в них сеточных функций. При выводе разностных уравнений предполагается, что потоки частиц, скорости дрейфа и коэффициенты диффузии постоянны
Рис. 1. Схема нерегулярной прямоугольной сетки
между узлами разностной сетки. Интегрирование дифференциальных выражений для потоков частиц (3) по контрольному объёму, выделенному на рис. 1 серым цветом, даёт их конечно-разностные аппроксимации, обеспечивающие устойчивость расчёта при больших градиентах напряжённости. При этом уравнения непрерывности для электронов и ионов принимают вид
?к+1 пк вг Т ^
• к+1 _ • к+1
А
Тй
А
3 к+1 + 3 к+1 3 к+1 _ 3 к+1 + А2 ^ + €г-1/2,з + Чд'+1/2 Чд-1/2 = в £ к
2&
Ас-
(5)
где
Д 6+1 — 6-1 д 0 + 1 — С?-1 АС = -о-, = "
2 ' ^ 2
Потоки частиц в радиальном и аксиальном направлениях вычислялись по экспоненциальным формулам
г к Рк+1 Б к ехР ( г к
^/¿+1/2,7 гг+1/2,5 р V
Зг.
к+1
€г + 1/2,5 А,
+ 1/2,5 Сг + 1/2
„к
г + 1/2,5
- рк+1 б?
¿ + 1,5' 1г + 1/2,5
ехр г
Ук
^¿ + 1/2,5
- 1
„к
Рк+1 £) к ехр ( ¿к ] _ Рк+1 Б к
'ч-; /г 5 + 1 /2 ^ V /г 5 + 1 /2 / "/¿5 + 1 /г !
3 г
к+1
4.5+1/2
,5+1/2 ^г,!+1 /2
к
,к+1 П к
г¿,j + l +1 /2
А<5 + 1/2
ехр г,
^ 1 г¿,j + 1 /2
1
Символы ¿к+1/2 и ,+1/2 определяются формулами
4+1/2,5 / ^к
^ + 1/2,5
(р?+1у - рУ, г
^,5 + 1/2
-5Ре/ ^^
б к
4,5+1/2
(р?3+1 - рУ
Верхний индекс к относится к моменту времени т? , при этом временной шаг состав ляет величину АТй = тк+1 — тк. Выражения для потоков з?+1
•к+1
; 3/ + в узлах сетки
^-1/2,5 Сг,5'-1/2
l
к
(г — 1/2, ]) и (г, ] — 1/2) имеют аналогичные представления. Входящие в них значения
потенциала , т^, 5^+1, коэффициентов переноса $г+1/2,3, $г,3+1/2, 0кг+1/2,3, А*3.+1/2, и источниковый член . в (5) берутся с предыдущего временного слоя . В силу локального выполнения законов сохранения для каждого контрольного объёма решение даже на грубой сетке удовлетворяет точным интегральным соотношениям. Полученная численная схема сохраняет положительные значения концентраций плазменных компонентов и электронной энергии.
Для вычисления электрического потенциала на к-м временном шаге использовалась конечно-разностная схема
а2 ^—^ ^—+ А (¿щ—^+7^-м, +
Д?» \ Д?г+1/2 Д?г—1/2 / 26 \ Д?г+1/2 Д?г—1/2
(%—^ — =Г(рк - 3 (6)
Д<л Д3+1/2 Дс.<-/2 / 1 7
Формулы для потока энергии электронов получаются аналогично формулам для потоков частиц в уравнениях непрерывности. Значения электронной температуры в узлах сетки находились из уравнения
з е^.1—-рк ■ /2 ■— ^11/2 ■
е» ,3 ' е» , 3_е» , У ег , 3 + а2 «»+1/2, 3_«г—1/2 ,3 +
5 Дтк Д?г
• к+1 + • к+1 + 1 _ + 1
,2^«г+1/2,3 + «г—1/2,3 , ^,3+1/2 ^,3—1/2 .
ч2 г+1/ 2,3 ■»»—1/2,3 ,
2£ + ДС3
+ ха2 Л'к+1 Ек+1 + ^+1 Ек+1 А +
2 Vе« г+1/2,3 «+1/2,3 •/е«г—1/2,3 1/2,3 /
+X 6е+1 + ^ек+1 ) + 0^. . = 0, (7)
2 V сг,3+1/2 с»,3+1/2 еСг,3 —1/2 43—1/2^ е»,3 ' 7
где компоненты потока энергии электронов в радиальном и аксиальном направлениях вычисляются по формулам
zki+1/2 j ^vk^+i^ exp (zj+1/2j) - «j D
3
^+1/2 j A«i+i/2 exp ( z*+i/2j) - 1
j*
zSj+1/2 k.j+i/2 exp (zjj+i/^ - tfjPkj+iD
Сг,3+1/2 Д?3+1/2 ехр (<3+1/2) — 1
Компоненты напряжённости электрического поля определяются как
тк+1 _ т к+1 т к+1 _ тк+1
Ек+1 = _ т^+1,7 Ек+1 = _
€г+1/2,3 Д , Сг,3+1/2 Д '
?г+1/2 ^3+1/2
Выражения для компонентов потока электронной энергии , и электри-
«г—1/2,3 с»,3—1/2
ческого поля Е|+1 ., Е^1 в узлах сетки (г — 1/2, ]) и (г, ] — 1/2) имеют аналогичные представления.
Для потоков электронов и ионов на границах расчётной области (на электродах и боковой стенке камеры) использовалась формула линейной интерполяции, позволяющая выразить поток на границе через значения потока внутри области, рассчитанные по экспоненциальной схеме. Например, для верхнего электрода соответствующая интерполяционная формула имеет вид
. 3 к+1 — 3 к+1
• к+1 = 3 к+1 , АС5+1/2 гC¿.J^+l/2 Згz¿.j+l = 3Ч:т/2 2 ' Ал. .
Это обеспечивает устойчивость расчёта потока на границе при больших градиентах электрического поля.
Начальные распределения концентраций электронов и ионов выбирались в виде аналитического решения, полученного в диффузионно-дрейфовом приближении [1]: — для ВЧ-разряда между плоскостями:
П
р0 = рР = 2sin (nZ )
— для ВЧ-разряда в цилиндрической области:
п
Р0 = Р0 = 2.3162 J (2.405£)sin(nC).
Начальное распределение электронной температуры полагалось постоянным по объёму реактора.
Численное моделирование 2D ВЧ-разряда даже на сравнительно грубых сетках требует больших объёмов вычислений и занимает значительное время. Для ускорения расчётов было выполнено распараллеливание вычислительного алгоритма на многопроцессорной системе с применением MPI. Для распараллеливания решения задачи использовалась декомпозиция расчётной области, что позволило однородно распределить данные для вычислений по процессорам. В [12] для случая трёхмерного уравнения конвективно-диффузионного переноса примеси от точечного источника, однотипного по структуре с уравнениями (1), было показано, что двумерная декомпозиция является более эффективной по сравнению с одномерной. На основании этого в предлагаемом алгоритме декомпозиция расчётной области проводилась по обеим координатам. Каждая из прямоугольных подобластей граничила не более чем с четырьмя соседними. Обмен между процессорами происходит только теми значениями j, р^1, Рр+/, , которые находятся на соприкасающихся границах подобластей, т. е. между минимальным количеством узлов сеточного шаблона, принадлежащих различным процессорам. В каждой подобласти рассматриваются четыре уравнения: для электрического потенциала, концентраций электронов и ионов, электронной температуры. Поэтому совместно с декомпозицией расчётной области использовалось распараллеливание по уравнениям (для отдельной подобласти на четыре процессора). Поскольку в численной схеме концентрации плазменных компонентов и электронная температура зависят от потенциала, то первоначально на одном процессоре в каждой подобласти решается конечно-разностное уравнение (6). Для перехода от момента времени Tk к моменту = Tk + ATk найденные значения потенциала j в отдельной подобласти передаются трём другим процессорам в этой же подобласти для решения конечно-разностных уравнений (5), (7). Значения Ре+/, Рр+/ и вычисляются в каждой подобласти одновременно на новом временном слое. Далее происходит обмен данными между четырьмя процессорами из одной
подобласти и данными между областями. Каждое из уравнений в сеточной области записывается в виде системы линейных алгебраических уравнений, имеющей матрицу большой размерности, но с разреженно-упорядоченной структурой. Для хранения таких матриц использовался метод типа CSR, позволяющий хранить несимметричные матрицы произвольной структуры [13]. В этом случае разреженная матрица хранится в виде одномерного массива, содержащего все ненулевые элементы, перечисленные в строчном порядке. Для решения таких СЛАУ сильная разреженность матрицы (сравнительно небольшое количество ненулевых элементов в строке или столбце) позволяет эффективно использовать классические итерационные и проекционные методы (с пред-обуславливанием) на подпространствах Крылова. В работе полученные СЛАУ решались итерационным методом Гаусса — Зейделя, показавшим хорошую эффективность.
В численных расчётах собственно плазмохимического травления применялась разработанная ранее авторами [14-16] математическая модель ПХР, в которой течение газовой смеси описывается уравнениями Навье — Стокса в приближении Буссинеска. Распределения концентраций частиц находятся из системы уравнений конвективно-диффузионного переноса. Пространственные распределения плотности электронов, определяющие скорости генерации химически активных частиц, брались из расчётов ВЧ-разряда в гидродинамическом приближении по описанному алгоритму. Плазмохимическая кинетика включала реакции диссоциации электронным ударом различных компонентов газовой смеси и реакции рекомбинации полученных атомов и радикалов с участием третьего тела. На поверхности образца рассматривались гетерогенные реакции адсорбции некоторых радикалов. Уравнения гидродинамики и конвективно-диффузионного переноса решались с использованием итерационной конечно-разностной схемы со стабилизирующей поправкой.
2. Результаты и обсуждение
В расчётах рассматривался плазмохимический реактор радиальной схемы при рабочем давлении p = 0.5 Торр и расходах газа Q = 200 ^ 800 см3/мин [14, 15].
Реакционная камера имела радиус электродов rt = 30 см и межэлектродное расстояние L = 3.5 см. Электроды предполагались термостатироваными при температуре Ts = 300 K. Состав рабочей смеси CF4/O2 изменялся в широких пределах. ВЧ-разряд инициировался на рабочей частоте f = 13.56 МГц при напряжении на электродах
= 110 В. Другие параметры были следующими: плотность газа n = 3 • 1016 см-3, средняя плотность электронов ne0 = 6 • 109 см-3, электронная температура на электродах и боковой стенке Tes = 0.5 эВ.
Характерные значения величин для обезразмеривания, взятые из работ [3, 4, 6], где рассматривались рабочие газы, близкие к исследуемому, и справочника [17], приведены ниже:
^e0 ^p0 De0 Dp0 k;0 E; h; he0
2 • 105 см2/В • с 2 • 103 см2/В • с 106 см2/с 102 см2/с 2.5 • 10-6 см3/с 24 эВ 15.578 эВ 15 эВ
Критерии и параметры, использованные в расчётах, были следующими:
T Pe Pp Da х Г вР Y 168.49 92 9200 1.6- 104 30.66 195.52 104 0.046
Расчёты ВЧ-разряда проводились на регулярной сетке 40 х 20 по пространственным переменным, использованной во всех вычислениях. Моделирование ВЧ-разряда с гармоническим током разряда требует решения нестационарных уравнений переноса зарядов и баланса энергии электронов. В расчётах рассматривалось ~ 104 периодов ВЧ-разряда, по истечении которых колебания электронной плотности и температуры выходили на периодический режим. Из-за сильной нелинейности исходных уравнений, включающих решение уравнения Пуассона, полученные неявные численные схемы являются условно устойчивыми. В этой связи шаг по времени АТк выбирался постоянным и достаточно малым — в диапазоне ~ 5 ■ 10-6 ^ 10-4.
Расчёт аксиально-симметричного ВЧ-разряда в двумерной постановке даже при малых удельных энерговкладах является трудоёмкой вычислительной задачей, решение которой занимает значительное время. При расчёте разряда в пространственной области применение параллельного вычислительного алгоритма, основанного на покомпонентном расщеплении системы уравнений и декомпозиции расчётной области, позволяет эффективно снизить время вычислений. В численных расчётах исходная область разбивалась на четыре равных подобласти, содержащие одинаковое количество узлов. Сравнение времени выполнения распараллеленной программы на 16 процессорах с временем выполнения исходной программы показало ускорение примерно в пять раз при эффективности распараллеливания 31 %. Эффективным ускоряющим фактором здесь является декомпозиция расчётной области. Теоретическая оценка ускорения при двумерной декомпозиции без учёта покомпонентного расщепления системы уравнений на квадратной сетке выражается формулой [12]
^Р/(1 + 4£ N),
где р — количество процессоров; ¿а, ¿0 — время выполнения арифметических операций и передачи данных между процессорами соответственно; N — размерность квадратной сетки ^ = N = N, в нашем случае N « 14). Поскольку решение СЛАУ занимает значительное время ¿а ^ ¿0, то с увеличением числа разбиений (задействованных процессоров р) ускорение растет практически линейно с увеличением р. Основным лимитирующим фактором, снижающим ускорение, является последовательность расчёта "потенциал — концентрации плазменных компонентов и электронная температура", которая заложена в численной схеме. Другой существенный лимитирующий фактор — ограничение шага по времени. С уменьшением шагов сетки А^.+х/2, А^5+1/2 (увеличением размерности N или использованием сгущающейся сетки для описания тонких приэлектродных слоёв) для обеспечения устойчивости счёта шаг по времени АТк также необходимо уменьшать. Снижение шага по времени пропорционально увеличивает общее время расчёта, что может перекрыть эффект от ускорения.
Тестирование программы проводилось на основе сравнения численного решения с формулой (8), полученной аналитически в диффузионно-дрейфовом приближении для ВЧ-разряда между плоскостями. Чтобы постановка задачи соответствовала тестовой, вместо краевого условия з'ес = —73РС на катоде (0 < £ < £ =1) бралось условие ре = 0. На рис. 2 пунктирной линией представлена зависимость электронной плотности от высоты на оси реактора, полученная в численных расчётах. Там же сплошной линией показана зависимость (8). Приведённые данные показывают качественное совпадение сопоставляемых результатов. При плоской конфигурации электродов расчётная электронная плотность однородна вдоль радиальной координаты, а её профиль
на оси симметрии практически совпадает с синусоидальной зависимостью от высоты реактора. Для количественной оценки различия электронных плотностей рассматривалось максимальное локальное отклонение. Найдено, что его значение не превышает величины 0.05, что находится в пределах ошибки аппроксимации.
На рис. 3 представлены распределения электронных плотностей в плазмохимиче-ском реакторе травления для плоской и цилиндрической конфигурации электродов. В обоих случаях электронная плотность максимальна в центре реактора и монотонно падает до нуля с приближением к электродам. При разряде в реакторе цилиндрической конфигурации электронная плотность снижается также и в радиальном направлении (см. рис. 3, б). Максимальная плотность электронов (в сечении г = 3.25 см) монотонно уменьшается вдоль радиуса на 0.3, 2.3, 17.1, 72% при г = 8.14, 14.57, 21.00, 27.43 см соответственно. Если не рассматривать крайнее значение г = 27.43 см, поскольку образцы располагаются на некотором расстоянии от внешней кромки электрода, снижение расчётного распределения электронной плотности по радиусу реактора составляет ~ 17 %, что должно существенно повлиять на однородность обработки образцов.
Отношение характерных времен для электронов и ионов составляет величину ~ 104, что соответствует существенно большей подвижности электронов по сравнению с иона-
Рис. 2. Зависимость электронной плотности от высоты на оси реактора: 1 — расчёт, 2 — 8т(п(")
a
_I__I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_i_L_I_
О 5 10 15 20 25 г, см
б
_I__I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_Ш_I_
0 5 10 15 20 25 г, см
Рис. 3. Изолинии электронной плотности пе ■ 10-10 см-3 в плазмохимическом реакторе: а — плоская, б — цилиндрическая конфигурации
ми. В расчётах это приводит к тому, что в высокочастотном электрическом поле движутся только электроны, а ионы остаются неподвижными. Сравнение расчётной электронной плотности в различные моменты периода колебаний напряжения показали, что концентрация электронов изменяется только вблизи электродов. Распределения ионных плотностей повторяют структуры распределений электронных плотностей и не реагируют на высокочастотное поле.
Модели физико-химической кинетики в зоне ВЧ-разряда включают значительное число реакций, инициируемых бомбардировкой электронами, в том числе ионизацию и диссоциацию молекул и радикалов, обменные реакции, реакции возбуждения. Основной трудностью при рассмотрении подробных кинетик является определение констант скоростей этих процессов. При известном распределении электронной температуры скорости кинетических процессов с участием электронов вычисляются по формуле
где е, о (е) — пороговая энергия и дифференциальное сечение соответствующего процесса, к — постоянная Больцмана. Расчёты показали, что электронная температура в объёме ПХР фактически постоянна, исключая узкие приэлектродные слои, где она резко меняется, но концентрация электронов мала. Среднее по периоду значение электронной температуры для обеих конфигураций ВЧ-разряда, которое принималось в последующих расчётах констант скоростей реакций, составляет ~ 5.9 эВ. Следует отметить, что эта величина, рассчитанная в рамках гидродинамической модели, оказалась близкой к полученной из решения кинетического уравнения Больцмана для электронов [18].
Для иллюстрации влияния электронной плотности на перенос нейтральных частиц, скорость и однородность травления образцов рассматривается случай, когда поток газовой смеси направлен к центру реактора. Расчётный режим характеризуется следующими параметрами: р = 0.5 Торр, Q = 200 см3/мин, Т = 300 К, содержание 02 в СГ4/02 составляет 20%. При других расчётных режимах результаты имеют аналогичный характер.
На рис. 4 приведены распределения концентрации фтора, соответствующие электронным плотностям, представленным на рис. 3. Несмотря на интенсивную конвекцию в радиальном плазмохимическом реакторе, основным механизмом переноса, определяющим процесс травления, является концентрационная диффузия активных частиц к поверхности обрабатываемых образцов. При распределении электронной плотности, показанной на рис. 3, а, концентрация фтора монотонно возрастает на входе в зону разряда и достигает максимума в диапазоне г =14 ^ 21 см, г = 3 ^ 5 см (рис. 4, а). Далее ввиду усиливающегося действия конвективного переноса концентрация фтора слабо понижается к центру реактора. При ВЧ-разряде в реакторе цилиндрической конфигурации (см. рис. 3, б) распределение концентрации фтора имеет подобный вид (см. рис. 4, б). Однако максимум концентрации фтора сдвинут ближе к центру реактора в зону координат г =11 ^ 16 см, г = 3 ^ 4.5 см, что связано с общим уменьшением электронной плотности и, как следствие, снижением производства активных частиц на периферии реактора.
На рис. 5 представлены зависимости скорости спонтанного травления от радиуса образца для двух конфигураций ВЧ-разряда и различного содержания 02 в СР4. При рассмотрении ВЧ-разряда между плоскостями скорость травления образца практиче-
о
ски не зависит от радиальной координаты (см. кривые 1-3). Однако из-за накопления активных частиц вблизи центра реактора скорость травления слабо повышается в диапазоне г = 14 ^ 21 см. Конкуренция между процессами образования и переноса фтора вблизи внешней кромки нижнего электрода понижает скорость травления образца в этой области. При ВЧ-разряде в реакторе цилиндрической конфигурации скорость травления почти постоянна в диапазоне г = 2 ^ 16 см и далее монотонно снижается в диапазоне г =16 ^ 27.5 см, что объясняется общим уменьшением электронной плотности на внешней кромке нижнего электрода и, как следствие, снижением скорости генерации активных частиц (см. кривые 4 - 6). Для обеих конфигураций ВЧ-разряда с увеличением содержания 02 с 10 до 30 % скорость травления у внешней кромки образца значительно понижается, а её профиль становится более неравномерным. Различие между зависимостями средних скоростей травления от содержания 02 в СР4 для двух конфигураций ВЧ-разряда не превышает 4.7%, что объясняется одинаковой в обоих случаях мощностью разряда.
Рис. 4. Изолинии концентраций фтора С • 109, моль/см3, в радиальном плазмохимическом реакторе: а — ВЧ-разряд между плоскостями, б — ВЧ-разряд в цилиндрическом объёме. Направление подачи смеси — к центру реактора
V,
А/мин 3000
2000
1500
1«=-Д- - - —1 ^-- - А _ -- -Д- -- ^ - —д
— — -^ X. N А
& - - < -© Ч> э - О
:
Л о £ II II мП
43- - -1111 О --Е II ]-- мм - В III -о 1111 1111 э -1111 II
- -О- - 1
- О ~ 2
- -¿^ - 3
-О- 4
-<3- 5
-Ь- б
10 12 14 16 18 20 22 24 г, см
Рис. 5. Распределение скорости травления у8 вдоль радиуса образца: 1, 2, 3 — ВЧ-разряд между плоскостями, 3, 4, 5 — ВЧ-разряд в реакторе цилиндрической конфигурации; 1, 4 — 10% 02; 2, 5 — 20% 02; 3, 6 — 30% 02. Другие параметры см. на рис. 3
Для количественного описания неоднородности травления вычислялся индекс неоднородности травления по формуле
т = —
где 1^1,1^2,^ — минимальное, максимальное и среднее значения спонтанной скорости травления. С увеличением содержания 02 в СР4 до 25 % в ВЧ-разряде между двумя поверхностями индекс неоднородности достигает 9.7%. Однако на практике обрабатываемые образцы размещаются на некотором расстоянии от внутренней и внешней кромок нижнего электрода. Заметим, что наибольшая неоднородность травления локализуется вблизи внешней кромки образца. Поэтому снижением скорости травления в этой зоне можно пренебречь, и неоднородность травления оказывается существенно ниже. При рассмотрении ВЧ-разряда в цилиндрической области индекс неоднородности увеличивается до 16.3 %, что объясняется снижением скорости травления в диапазоне г = 17-27.5 см. "Плоская" конфигурация ВЧ-разряда обеспечивает более равномерное производство активных частиц и лучшую однородность травления.
Заключение
Разработан алгоритм численного моделирования ВЧ-разряда в гидродинамическом приближении. Для решения уравнений непрерывности электронов и ионов и уравнения энергии электронов использована обобщённая на двумерный случай неявная экспоненциальная разностная схема, которая обеспечивает положительную определённость расчётных значений электронной температуры и концентраций плазменных компонентов. На основе этой схемы рассчитаны характеристики разряда в плазмохимическом реакторе травления радиальной схемы для плоскопараллельной и цилиндрической конфигураций ВЧ-электродов. Опыт использования предложенного алгоритма для гидродинамической модели ВЧ-разряда в расчёте режимов травления в ПХР показал, что время отдельного расчёта параллельной программой составляет 3 - 5 сут. Поэтому из-за большой вычислительной трудоёмкости численное моделирование ВЧ-разряда следует рассматривать как отдельное самостоятельное исследование, а в оптимизационных расчётах ПХР для электронной плотности и температуры целесообразно применять аналитические модели, учитывающие результаты численного моделирования или экспериментальные данные.
Список литературы
[1] Cherrington B.E. Gaseous Electronics and Gas Lasers. Oxford: Pergamon Press, 1980.
[2] Aydil E.E., Eoonomou D.J. Modeling of plasma etching reactors including wafer heating effects // J. of the Electrochem. Soc. 1993. Vol. 140. P. 1471-1481.
[3] Graves D.B., Jensen K.F. A continuum model of DC and RF discharges // IEEE Trans. on Plasma Sci. 1986. Vol. PS-14. P. 78-91.
[4] Lymberopoulos D.P., Economou D.J. Fluid simulation of glow discharges: Effect of metastable atoms in argon // J. of Appl. Phys. 1993. Vol. 73. P. 3668-3679.
[5] Глдияк Г.В., Швейгерт В.А., УуэмАА О.У. Математическое моделирование тлеющего газового разряда // Изв. СО АН СССР. Серия техн. наук. 1988. № 21, вып. 6. С. 41-47.
[6] Швейгерт В.А. Высокочастотный разряд низкого давления в электроотрицательных газах. Новосибирск, 1990 (Препр. РАН. Сиб. отд-ние. ИТПМ. № 8-90).
[7] Тихонов А.Н., Самарский А.А. Об одной наилучшей однородной разностной схеме // Докл. АН СССР. 1959. Т. 124. С. 779-782.
[8] Дулан Э., Миллер Дж., Шилдерс У. Равномерные численные методы решения задач с пограничным слоем. М.: Мир, 1983.
[9] Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984.
[10] Scharfetter D.L., Gummel H.K. Large-signal analysis of a silicon read diode oscillator // IEEE Trans. on Electron Devices. 1969. Vol. ED-16. P. 64-77.
[11] Ting Wei Tang. Extension of the Scharfetter-Gummel algoritm to the energy balance equation // Ibid. 1984. Vol. ED-31. P. 1912-1914.
[12] Старченко А.В., Ибраев Г.М. Опыт создания вычислительного кластера на базе кластерных систем Томского научного центра // Четвёртая Сибирская школа-семинар по параллельным и высокопроизводительным вычислениям. Томск: Дельтаплан, 2008. С. 61-78.
[13] Писсанецки С. Технология разреженных матриц. М.: Мир, 1988.
[14] Григорьев Ю.Н., Горобчук А.Г. Эффекты неизотермичности в плазмохимическом реакторе травления // Микроэлектроника. 1998. Т. 27, № 4. С. 294-303.
[15] Григорьев Ю.Н., Горобчук А.Г. Особенности интенсификации травления кремния в CF4/O2 // Там же. 2007. Т. 36, № 4. С. 283-294.
[16] Grigoryev Yu.N., Gorobchuk A.G. Numerical simulation of plasma-chemical processing semiconductors // Micro Electronic and Mechanical Systems / Ed. by Kenichi Takahata. In-Tech Education and Publ., 2009. P. 185-210.
[17] Физические величины. Справочник / Под ред. И.С. Григорьева, Е.З. Мейлихова. М.: Энергоатомиздат, 1991.
[18] Meeks E., Vosen S.R., Shon J.W. et al. Results From Modeling and Simulation of Chemical Downstream Etch Systems. Sandia Rep. SAND96-8241 UC-401. 1996.
Поступила в 'редакцию 6 мая 2013 г.
