Физика плазмы, 2019, T. 45, № 8, стр. 755-768

Моделирование индукционного разряда в аргоне с учетом протока и неоднородности температуры газа

А. Н. Кропоткин a, Д. Г. Волошин a*

a МГУ им. М.В. Ломоносова, НИИ ядерной физики им. Д.В. Скобельцына, (НИИЯФ МГУ)
Москва, Россия

* E-mail: dvoloshin@mics.msu.su

Поступила в редакцию 19.12.2018
После доработки 30.01.2019
Принята к публикации 07.02.2019

Полный текст (PDF)

Аннотация

С целью развития методов контроля параметров низкотемпературной неравновесной плазмы разработана численная модель индукционного высокочастотного разряда в аргоне в диффузионно-дрейфовом приближении и проведено теоретическое исследование транспорта ионов на поверхность обрабатываемого материала. Построенная модель была протестирована на имеющихся экспериментальных и теоретических данных. Расчеты проводились для индукционного разряда в аргоне с характерными для современных реакторов параметрами: частота 13.56 МГц при давлении в камере 10 мТорр. Модель использовалась для определения плотности плазмы, температуры электронов и потока ионов на обрабатываемый материал. Получены данные о температуре газа в зависимости от вложенной мощности и о параметрах разряда в зависимости от протока газа.

1. ВВЕДЕНИЕ

Индукционные плазменные реакторы широко используются в современной технологии производства микроэлектроники [1], в медицинских устройствах [2], как ионные источники для ускорителей [3] и термоядерного синтеза [4]. Особенности индукционного разряда – возможность получения высоких плотностей плазмы (1011–1012 см–3) даже при низких давлениях газа (около 10 мТорр), и возможность работы без прямого контакта электродов с плазмой. Применение индукционных разрядов в микроэлектронике включает в себя процессы осаждения и травления полупроводников и металлов [5]. В настоящий момент идет исследование применения индукционных разрядов в новых технологических процессах атомно-слоевого травления [6], возможности работы в импульсном режиме [7] и исследования элементарных процессов в подобных установках [8]. Вследствие этого возрос интерес к построению моделей таких реакторов [9, 10].

Актуальной темой в микроэлектронике является изучение материалов с низкой диэлектрической проницаемостью (low-k) [1113]. Использование таких материалов в качестве межслойных диэлектриков в микросхемах представляет собой один из подходов для повышения объемной плотности элементов в микроэлектронных устройствах. Также эти материалы позволяют уменьшить емкость между слоями в микрочипах и тем самым снизить задержку прохождения сигнала. Изучение воздействия плазмы на low-k материалы является актуальной задачей и требует прецизионного контроля параметров плазмы: потоков и спектров ионов, потоков нейтральных частиц. Данные об этих параметрах плазмы могут быть получены при условии использования адекватных численных моделей для описания процессов в плазме индукционных и емкостных высокочастотных (ВЧ) разрядов.

Разработка самосогласованных численных моделей газовых разрядов – комплексная задача. Эти модели включают в себя решение кинетических уравнений для нейтральных и заряженных частиц, электромагнитные уравнения Максвелла, уравнение Навье–Стокса, если необходимо учесть проток газа, уравнения теплового баланса для нейтральных частиц, чтобы учесть нагрев газа, а также большое число объемных и поверхностных химических реакций. Моделирование разрядов активно используется для изучения фундаментальных задач физики плазмы и помогает в создании новых плазменных установок, в том числе и промышленных.

Кинетическое моделирование на основе метода “частиц в ячейке” позволяет учитывать нелокальную кинетику частиц без дополнительных упрощающих предположений, однако требует значительных вычислительных ресурсов. Поэтому в случае моделирования индукционных разрядов чистое кинетическое моделирование проводилось для режимов с низкой плотностью плазмы [14], или и с использованием сложных численных моделей [15], позволяющих ускорить расчет для больших плотностей плазмы. В большинстве случаев в зависимости от задач выбирается приближение, которое помогает упростить модель [1619] и снизить объем необходимых вычислений. Подходы активно развиваются на протяжении последних десятилетий и на данный момент ведутся активные исследования низкотемпературной плазмы с помощью моделирования [9].

В данной работе рассматривается построение самосогласованной гидродинамической модели индукционного разряда. В данном подходе из кинетического уравнения Больцмана для заряженных частиц получают систему гидродинамических уравнений, или так называемые моменты функции распределения электронов по энергии: уравнения непрерывности для плотности, сохранения импульса и сохранения энергии [2025]. Уравнения непрерывности также записываются и для нейтральных частиц. Далее эти уравнения решаются самосогласованно с уравнениями Максвелла, теплового баланса и Навье–Стокса для потока газа. В зависимости от условий разряда может быть сделано предположение о Максвелловском виде функции распределения электронов по энергии (ФРЭЭ) или вид ФРЭЭ может находиться из уравнения Больцмана в двучленном приближении для заданного значения электрического поля [26, 27]. Недостатком такой гидродинамической модели является отсутствие данных о функции распределения ионов по энергии (ФРИЭ). В данной работе для преодоления этого недостатка модель индукционного разряда была дополнена кинетической моделью для расчета ФРИЭ в заданном электрическом поле из гидродинамической модели с учетом столкновений с нейтральными частицами методом Монте-Карло.

Таким образом представленная в данной работе модель отличается от многих моделей индукционного разряда, которые не включают в себя отдельных уравнений для температуры и протока газа, как, например, [24, 28]. А также от моделей с учетом неоднородности температуры и протока газа, но без возможности кинетического расчета ФРИЭ [22].

Экспериментальные зондовые исследования показывают, что во многих случаях ФРЭЭ в индукционных разрядах близка по форме к Максвелловской функции [2931]. Это позволяет получать точные результаты при использовании гидродинамических моделей в предположении Максвелловской ФРЭЭ [22, 24, 28], которые требуют значительно меньше времени для достижения стационарного решения, чем кинетические подходы, особенно когда моделируется реактор, геометрию которого нельзя свести к одномерной (что типично для индукционных реакторов).

Статья организована следующим образом. В разд. 2 описывается численная модель, состоящая из гидродинамической модели плазмы с учетом температуры и протока газа, а также кинетической модели движения ионов в заданном поле с учетом столкновений для определения энергетического спектра ионов на нижнем электроде. В разд. 3 рассматриваются и обсуждаются основные результаты расчетов. Приводятся данные о тестировании гидродинамической модели с помощью численных и экспериментальных работ других авторов. Представлен результат расчета разряда в условиях эксперимента [11] и сравнение с экспериментально измеренными параметрами. Также в разд. 3 приведены результаты по зависимости температуры газа от вложенной в разряд мощности и зависимости параметров разряда от величины расхода газа для геометрии разряда [11]. В разд. 4 даются выводы работы.

2. ЧИСЛЕННАЯ МОДЕЛЬ

2.1. Гидродинамическая модель плазмы индукционного разряда с учетом температуры и протока газа

Плотность плазмы, потоки ионов на электрод, и температура электронов находятся из решения системы гидродинамических уравнений. В эту систему входят следующие уравнения.

Уравнение непрерывности для плотностей частиц (электронов, ионов, нейтральных частиц)

(1)
$\frac{{\partial n}}{{\partial t}} + \nabla \cdot \left( {n{\mathbf{v}}} \right) = \mathop \sum \limits_j {{R}_{j}},$
где n – плотность и v – направленная скорость частиц. В правой части стоит сумма скоростей рождения и уничтожения частиц в результате j-й реакции.

Закон сохранения импульса частиц

(2)
$\frac{\partial }{{\partial t}}\left( {nm{\mathbf{v}}} \right) + \nabla \cdot \left( {nm{\mathbf{vv}}} \right) = - \nabla P + n{\mathbf{F}} - nm{\mathbf{v}}{{\nu }_{m}}$
можно записать в виде
(3)
$mn\left( {\frac{\partial }{{\partial t}}{\mathbf{v}} + ({\mathbf{v}} \cdot \nabla ){\mathbf{v}}} \right) = - \nabla P + n{\mathbf{F}} - nm{\mathbf{v}}{{\nu }_{m}},$
где νm – транспортная частота, P – парциальное давление рассматриваемой компоненты, F – внешняя сила. Если распределение частиц по энергиям можно считать близким к распределению Максвелла, то давление в таком случае полагают изотропным, P = nκT, где κ – постоянная Больцмана, а T – температура частиц. В большинстве случаев членами в левой части уравнения (2) для электронов можно пренебречь [21]. Первый член пренебрежимо мал, если v плавно меняется во времени. Вторым, инерционным, в случае электронов можно пренебречь, если он много меньше члена, отвечающего за столкновения [32], что верно при давлениях выше порядка 5 мТорр.

Ионы, однако, намного массивней и не могут также быстро реагировать на изменения электрического поля. В таком случае для ионов необходимо решать полное уравнение закона сохранения импульса или вводить эффективное поле [33]. Проведенные оценки показали, что в случае отсутствия дополнительного ВЧ-напряжения на электроде в нижней камере в конфигурации разряда [11], рассматриваемого в разделе 3.2, учет полного уравнения для момента импульса ионов незначительно влияет на результаты расчетов. Эффект учета полного уравнения для момента импульса ионов будет рассмотрен в дальнейшей работе по моделированию индукционного разряда с дополнительным ВЧ-напряжением на электроде.

Таким образом, в данной модели для электронов и ионов используется диффузионно-дрейфовое приближение [25]

(4)
${\mathbf{J}} = \mu n{\mathbf{E}} - D\nabla n,$
где $\mu = q{\text{/}}(m{{\nu }_{m}})$ – подвижность частиц, q – заряд частицы, а $D = \kappa T{\text{/}}(m{{\nu }_{m}})$ – коэффициент диффузии.

Граничное условие для потоков электронов и ионов [21]

(5)
${\mathbf{n}} \cdot {{{\mathbf{J}}}_{{\mathbf{e}}}} = \frac{1}{2}{{v}_{{th}}}{{n}_{e}},\quad {\mathbf{n}} \cdot {{{\mathbf{J}}}_{i}} = {{v}_{{bohm}}}{{n}_{i}},$
где ${{v}_{{th}}}$ – тепловая скорость электронов, ${{v}_{{bohm}}} = \sqrt {e{{T}_{e}}{\text{/}}{{M}_{i}}} $ – Бомовская скорость ионов.

Закон сохранения энергии для электронов в общем виде

(6)
$\begin{gathered} \frac{\partial }{{\partial t}}\left( {n\varepsilon } \right) + \nabla \cdot \left( {n{\mathbf{v}}\varepsilon } \right) = - \nabla P \cdot {\mathbf{v}} + nm{\mathbf{h}} \cdot {\mathbf{v}} + \\ \, + \nabla K \cdot \nabla T - \mathop \sum \limits_j {{R}_{j}}{{H}_{j}}, \\ \end{gathered} $
где K – коэффициент теплопроводности и Hj – энтальпия j-й реакции. Учитывая, что для распределения Максвелла ε = 3κT/2, а также предположения, сделанные при выводе диффузионно-дрейфового приближения, получим уравнение для кинетической энергии хаотического движения частиц
(7)
$\frac{\partial }{{\partial t}}\left( {\frac{3}{2}n\kappa T} \right) + \nabla \cdot {\mathbf{q}} + e{{{\mathbf{J}}}_{e}} \cdot {\mathbf{E}} + \mathop \sum \limits_j {{R}_{j}}{{H}_{j}} = 0,$
(8)
${\mathbf{q}} = - K\nabla T + \frac{5}{2}\kappa T{{{\mathbf{J}}}_{e}},$
где q тепловой поток.

Граничное условие для концентрации энергии электронов [21]

(9)
${\mathbf{n}} \cdot {\mathbf{q}} = \frac{5}{6}{{v}_{{th}}}{{n}_{\varepsilon }}.$

Для ионов уравнение для кинетической энергии хаотического движения не рассчитывалось. Вместо этого в соответствии с [34] было сделано предположение, что температура ионов постоянна во всей камере и равняется порядка 0.1 эВ.

Константы электронных реакций рассчитывались по формуле

(10)
${{k}_{j}} = \mathop \smallint \limits_0^\infty {{\sigma }_{j}}\left( \varepsilon \right){{u}_{e}}\left( \varepsilon \right)f(\varepsilon )d\varepsilon ,$
где σj – это сечение рассеяния для j-го процесса, f(ε) – функция распределения электронов по энергиям (ФРЭЭ), ue – модуль скорости электронов, а $\varepsilon = {{m}_{e}}u_{e}^{2}{\text{/}}2$ – энергия электронов. В табл. 1 приведена схема реакций, используемых при расчетах. ФРЭЭ нормировалась на единицу. Существует два основных подхода для расчета констант реакций (10). В первом ФРЭЭ полагают равной распределению Максвелла или распределению Дрювестейна и по известному набору сечений рассчитывают kj как функции средней энергии. Во втором подходе используют отдельный решатель уравнения Больцмана в постоянном поле в двучленном приближении для расчета ФРЭЭ, и затем находят kj как функцию средней энергии ε.

Таблица 1.

Схема реакций. Электронные возбуждения и нейтральная химия

# Реакция Процесс Порог энергии, скорость процесса Ссылка
(R1) ${\text{Ar}} + e \to {\text{Ar}} + e$ Упругое рассеяние [35]
(R2) ${\text{Ar}} + e \leftrightarrow {\text{Ar*}} + e$ Возбуждение 11.56 [35]
(R3) ${\text{Ar}} + e \to {\text{A}}{{{\text{r}}}^{ + }} + 2e$ Ионизация 15.8 [35]
(R4) ${\text{Ar*}} + e \to {\text{A}}{{{\text{r}}}^{ + }} + 2e$ Ступенчатая ионизация 4.14 [36]
(R5) ${\text{Ar*}} + {\text{Ar*}} \to {\text{A}}{{{\text{r}}}^{ + }} + e + {\text{Ar}}$ Хим. ионизация 6.2 × 10–10 [37]
(R6) ${\text{A}}{{{\text{r}}}^{ + }} \to {\text{Ar}}$ Гибель на стенке γ = 1  
(R7) ${\text{Ar}}* \to {\text{Ar}}$ Гибель на стенке γ = 0.8 [38]

Здесь введены обозначения: Ar – основное состояние аргона, Ar* – суммарное возбужденное состояние аргона (порог реакции 11.5 эВ), Ar+ – ион аргона (порог ионизации 15.8 эВ).

В данной работе ФРЭЭ предполагается максвелловской. Результаты экспериментальных измерений [30, 31] показывают справедливость такового предположения в исследуемом диапазоне давлений для плазмы Ar.

Распределение электромагнитных полей находится из решения уравнения Максвелла. В данной работе используется двумерная аксиально-симметричная геометрия, и в таких условиях уравнения Максвелла можно свести к уравнению [39]

(11)
$\frac{1}{{{{\mu }_{0}}{{\mu }_{r}}}}{{\nabla }^{2}}{\mathbf{A}} + {{\varepsilon }_{0}}{{\varepsilon }_{r}}{{\omega }^{2}}{\mathbf{A}} = - {{{\mathbf{J}}}_{s}},$
где A – магнитный потенциал, а Js – плотность тока, создаваемого внешним источником.

В качестве входного параметра в данной работе задавалась мощность, вкладываемая в катушки, ${{P}_{{coil}}} = {{V}_{{coil}}}I_{{coil}}^{*}{\text{/}}2$. Отсюда рассчитывается потенциал катушек, ${{I}_{{coil}}} = \int {{{{\mathbf{J}}}_{i}} \cdot \varphi dS} $ ток в катушках, где Ji – плотность тока, индуцируемого в плазме. Этот параметр определяет член Js в уравнении (11)

(12)
${{{\mathbf{J}}}_{s}} = \frac{{{{\sigma }_{p}}}}{{2\pi r}}{{V}_{{coil}}}{\mathbf{\varphi }},$
σp – комплексная проводимость плазмы (в случае максвелловской ФРЭЭ),
(13)
${{\sigma }_{p}} = \frac{{{{n}_{e}}{{e}^{2}}}}{{{{m}_{e}}({{\nu }_{m}} + i{{\omega }_{{RF}}})}}\quad,$
где me – масса электрона, νm – эффективная частота упругих соударений.

Все железные стенки камеры предполагаются заземленными, кроме плавающего электрода, и кварцевых стенок где

(14)
$ - {\mathbf{n}} \cdot {\mathbf{D}} = {{\rho }_{s}},\quad \frac{{\partial {{\rho }_{s}}}}{{\partial t}} = {\mathbf{n}} \cdot {{{\mathbf{J}}}_{i}} + {\mathbf{n}} \cdot {{{\mathbf{J}}}_{e}},$
где Ji и Je – потоки ионов и электронов соответственно. Граничное условие для векторного потенциала на внешних границах моделируемого объема
(15)
${\mathbf{n}} \times {\mathbf{A}} = 0,$
что означает магнитную изоляцию.

Для лучшего описания параметров разряда необходимо рассчитывать температуру нейтрального газа исходя из уравнения теплового баланса

(16)
$\rho {{C}_{p}}\frac{{\partial T}}{{\partial t}} + \rho {{C}_{p}}{\mathbf{u}} \cdot \nabla T + \nabla \cdot {\mathbf{q}} = Q,$
где ρ – плотность, Cp – теплоемкость при постоянном давлении, u – скорость, κ – коэффициент теплопроводности, Q – источники тепла, q – поток тепла, равный

(17)
${\mathbf{q}} = - \kappa \nabla T.$

Источник тепла от химических реакций $Q = $ $ = \sum\nolimits_j {{{Q}_{j}}} $, ${{Q}_{j}} = --{{H}_{j}}{{r}_{j}}$. Здесь rj – скорость реакции, ${{H}_{j}} = \sum\nolimits_i {{{v}_{{ij}}}{{h}_{i}}} $ – энтальпия реакции j, где ${{v}_{{ij}}}$ – стехиометрический коэффициент, а энтальпия hi записывается через термохимические полиномиальные коэффициенты [40]

(18)
${{h}_{i}} = R\left( {{{a}_{1}}T + \frac{{{{a}_{2}}{{T}^{2}}}}{2} + \frac{{{{a}_{3}}{{T}^{3}}}}{3} + \frac{{{{a}_{4}}{{T}^{4}}}}{4} + \frac{{{{a}_{5}}{{T}^{5}}}}{5} + {{a}_{6}}} \right),$
где R – универсальная газовая постоянная. Индекс i означает, что коэффициенты для возбужденных компонент или ионов могут отличаться от коэффициентов основного состояния. В этой работе данные для полиномиальных коэффициентов брались из базы данных NASA [40].

Граничные условия для температуры могут быть заданы в виде фиксированной температуры стенки, равенства нулю теплового потока (теплоизоляция стенки), либо в виде заданного потока

(19)
$ - {\mathbf{n}} \cdot {\mathbf{q}} = h({{T}_{{ext}}} - T),$
где h – коэффициент теплопередачи, Text – температура с внешней стороны стенки, при этом предполагается заданная внешняя температура (например, комнатная) на определенном расстоянии от стенок камеры. В данной работе это расстояние полагалось равным 0.1 м.

В промышленных и экспериментальных установках, в реакторе присутствует проток газа. Так как на поверхности обрабатываемого материала и в объеме реактора идут химические реакции, то обычно обеспечивают постоянный проток газа для удаления продуктов реакций и для поддержания концентрации реагентов. Для учета протока газа в модель включены уравнения Навье–Стокса. В стационарном случае и в предположении отсутствия турбулентности и объемных сил они имеют вид

(20)
$\begin{gathered} \rho \left( {{\mathbf{u}} \cdot \nabla } \right){\mathbf{u}} = - \nabla p + \\ \, + \nabla \left( {\mu \left( {\nabla \cdot {\mathbf{u}} + {{{\left( {\nabla \cdot {\mathbf{u}}} \right)}}^{T}}} \right) - \frac{2}{3}\mu \left( {\nabla \cdot {\mathbf{u}}} \right){\mathbf{I}}} \right), \\ \end{gathered} $
(20a)
$\rho \nabla \cdot {\mathbf{u}} = 0,$
где ρ – плотность газа, p – давление, µ – коэффициент динамической вязкости, верхний индекс T в правой части обозначает транспонирование, I – единичная матрица, ${\mathbf{u}} = 0$ – граничное условие для скорости нейтрального газа на стенках камеры.

Полученная в разделе система дифференциальных уравнений решалась методом конечных элементов с помощью пакета COMSOL [41].

2.2. Кинетический расчет функции распределения ионов по энергии

При описании взаимодействия плазмы с поверхностью важным параметром является распределение приходящих на поверхность частиц по энергии. Для расчета ФРИЭ гидродинамическая модель разряда была дополнена кинетической моделью движения ионов в заданном электрическом поле с учетом столкновений с нейтральными частицами методом Монте-Карло.

Заданное электрическое поле бралось из стационарного решения диффузионно-дрейфовой модели разряда. Динамика ионов рассчитывалась с учетом процессов упругого рассеяния и резонансной перезарядки на атомах Ar. В модели использованы сечения столкновений ионов Ar+ с атомами Ar из работы [42]. В процессе моделирования отслеживались траектории движения 105 частиц. После столкновения иона с нижним электродом его энергия записывалась в общий спектр.

3. РЕЗУЛЬТАТЫ РАСЧЕТОВ И ОБСУЖДЕНИЕ

3.1. Тестирование модели

Тестирование модели проводилось на основе имеющихся теоретических и экспериментальных данных для разряда в аргоне для разных геометрий установок.

Было проведено моделирование индукционного разряда в простой цилиндрической геометрии при следующих условиях: частота 6.78 МГц, давление газа 10 мТорр, мощность, вкладываемая в катушки 100 Вт, температура газа 500 К. Цилиндрическая камера разряда имела радиус 10 см и высоту 10.75 см. Антенна представляла собой плоскую катушку из пяти витков отделенную от камеры кварцевым стеклом толщиной 1.27 см. После проведения расчета были получены результаты, хорошо согласующиеся как с экспериментальными измерениями в данной геометрии [43], так и результатами моделирования из работы [24]. Сравнение основных характеристик разряда – аксиальных профилей плотности плазмы и температуры электронов, полученных в данной модели с экспериментальными результатами [43], представлено на рис. 1. Согласие расчетных и экспериментальных данных одновременно по плотности плазмы и температуре электронов говорит о возможности использования максвелловской формы ФРЭЭ в модели. При немаксвелловской форме ФРЭЭ необходимую плотность плазмы возможно получить только с отличной от экспериментальной температурой электронов.

Рис. 1.

Сравнение аксиальных распределений плотности плазмы и температуры электронов в расчете и эксперименте.

Также был проведен расчет для индукционного разряда в геометрии GEC (Gaseous Electronics Conference) [44]. Данная геометрия представляет собой стандартизированную систему для сравнения экспериментальных измерений и численных расчетов емкостных и индукционных разрядов, сделанных разными методами и разными научными группами в одинаковых условиях. Расчет индукционного разряда проводился для следующих параметров: разряд в аргоне на частоте 13.56 МГц при давлении 10 мТорр и температуре газа 500 К. Вложенная в разряд мощность составляла 77, 159 и 245 Вт в соответствии со значениями из [44]. Сравнение плотности N и потенциала плазмы ${{V}_{{pl}}}$ в центре разряда в зависимости от мощности представлено в табл. 2. Видно хорошее согласие расчетных данных с экспериментом.

Таблица 2.

Параметры плазмы для GEC-разряда

P, Вт N, 1011 см–3 эксперимент N, 1011 см–3 расчет ${{V}_{{pl}}}$, В эксперимент ${{V}_{{pl}}}$, В расчет
77 3.2 2.6 20 19.6
159 5.7 5.5 20.5 18.8
245 7.6 7.8 20.6 18.4

Таким образом, построенная модель была успешно протестирована на имеющихся экспериментальных и расчетных данных.

3.2. Расчет индукционного разряда вниз по потоку

В данном разделе рассмотрим моделирование разряда в экспериментальной установке из работы [11]. Данная установка представляет собой плазменный реактор, состоящий из двух частей: верхней камеры, в которой создается плазма индукционного разряда, и нижней камеры с дополнительным электродом c плазмой вниз по потоку из верхней камеры (“down stream plasma”). Схема установки приведена на рис. 2. Верхняя камера представляет собой кварцевый цилиндр диаметром 8 см и длиной 25 см. Снаружи посередине цилиндра размещается катушка из четырех витков медной проволоки. Плазма возбуждается на частоте 13.56 МГц. Нижняя камера представляет собой цилиндр из нержавеющей стали 35 см в диаметре и 10 см в высоту. Газ входит в камеру сверху кварцевой трубы и выходит снизу нижней камеры. Такая двухкамерная конфигурация разряда позволяет, при необходимости, отсекать плазму с помощью сетки и проводить исследования влияния как отдельных компонент из плазмы (радикалов, потока ВУФ (вакуумный ультрафиолет) фотонов), так и их совместного воздействия с ионами на образцы, располагаемые на электроде в нижней камере. Данная установка успешно использована для изучения механизмов воздействия плазмы на материалы с низкой диэлектрической проницаемостью [11, 45]. Однако для дальнейшего развития моделей взаимодействия плазмы с различными материалами необходимо описание состава и потоков частиц, прилетающих из плазмы, в том числе в сложных газовых смесях.

Рис. 2.

Схематичное изображение установки.

В данной работе полученные расчетные результаты в аргоне сравнивались с экспериментальными данными, приведенными в работе [11]. В эксперименте плотность плазмы измерялась зондом Hairpin, и поток ионов измерялся с помощью анализатора энергии ионов (расположение показано на рис. 2). Согласно условиям, приведенным в работе [11], расчет проводился в Ar на частоте 13.56 МГц, при давлении в камере 10 мТорр и вкладываемой мощности в катушки 100 Вт. Схема разбиения моделируемого объема на ячейки для расчета методом конечных элементов приведена на рис. 3.

Рис. 3.

Схема разбиения моделируемого объема на ячейки.

Пространственные распределения плотности плазмы, температуры электронов и потенциала электрического поля, полученные при моделировании, показаны на рис. 4–6 соответственно.

Рис. 4.

Распределение плотности плазмы: Ar, 13.56 МГц, 10 мТорр, 100 Вт. Верхняя и нижняя камеры имеют раздельную шкалу.

Рис. 5.

Распределение температуры электронов. Ar, 13.56 МГц, 10 мТорр, 100 Вт.

Рис. 6.

Распределение потенциала электрического поля. Ar, 13.56 МГц, 10 мТорр, 100 Вт.

Из распределения плотности плазмы видно, что область плотной плазмы располагается непосредственно около катушек с максимумом в районе $3.7 \times {{10}^{{11}}}$ см−3 и падает на два порядка в нижней камере до уровня порядка 109 см−3 около нижнего электрода. Вследствие такого градиента плазмы на рис. 4 верхняя и нижние камеры имеют отдельные обозначения (шкалы) для большей наглядности результатов. Также далее в разд. 3.3 на рис. 12 показан профиль плотности плазмы на оси в логарифмическом масштабе.

Температура электронов имеет максимум порядка 3.3 эВ в верхней камере в районе катушки и спадает до величины порядка 2 эВ в нижней камере.

Температура газа рассчитывалась исходя из уравнения теплового баланса (16), ее распределение приведено на рис. 7. Распределение газовой температуры сильно зависит от граничных условий на стенках камеры. Подробное обсуждение возможных вариантов граничных условий приведено в следующем разделе. На рис. 7 предполагается заданная температура стенок в нижней камере ${{T}_{{wall}}} = 405$ К, соответствующая экспериментальным наблюдениям, и поток тепла (19) в верхней камере.

Рис. 7.

Распределение температуры нейтрального газа.

В табл. 3 приведены результаты по плотности плазмы в нижней камере в точке соответствующей расположению датчика Hairpin на рис. 2 ($z = - 0.2$, $r = 0$), потоку ионов на нижний электрод в эксперименте и в различных вариантах расчета: а) без учета протока и с фиксированным значением температуры газа во всей камере, б) без учета протока с расчетом распределения температуры газа ${{T}_{{gas}}}$, в) с учетом протока и фиксированной температурой ${{T}_{{gas}}} = 500$ К, г) с учетом протока нейтрального газа и расчетом ${{T}_{{gas}}}$ (полная модель).

Таблица 3.

Сравнение значений плотности электронов и потока ионов с экспериментальными данными из работы [11]

  а) Без протока ${{T}_{{gas}}} = 500$ К б) Без протока Расчет ${{T}_{{gas}}}$ в) С протоком ${{T}_{{gas}}} = 500$ К г) С протоком Расчет ${{T}_{{gas}}}$ Эксперимент
ne, см–3 9.53 × 108 3.3 × 109 4.63 × 109 2.18 × 109 1.54 × 109
Fi, см–3 c–1 1.18 × 1014 3.86 × 1014 5.64 × 1013 2.36 × 1014 1.78 × 1014

Массовый расход газа был задан в соответствии с экспериментальным значением в работе [11] – 3 sccm (млн/мин). На рис. 8 приведено распределение скорости нейтрального газа.

Рис. 8.

Распределение скорости при массовом расходе газа 3 sccm.

Из приведенных в табл. 2 значений видно, что результаты расчетов в полной модели и экспериментальные данные согласуются (с учетом экспериментальной погрешности 30–50% для измерения плотности плазмы порядка 109 с помощью Hairpin зонда). При этом величина плотности плазмы и поток ионов в нижней камере зависят от используемого приближения для температуры и протока газа. При одинаковом подходе к расчету Tgas учет протока уменьшает значение плотности плазмы в нижней камере. Зависимость параметров плазмы от протока газа обсуждается в следующем разделе. Расчет температуры газа по сравнению с приближением постоянной температуры с учетом разогрева Tgas = 500 К увеличивает значение плотности плазмы в нижней камере. Это связано как с более плавными градиентами плотности плазмы в случае расчета температуры, так и с увеличением роли процессов ионизации (R4), (R5). В случае сценариев расчетов б), г) их вклад в суммарную скорость ионизации в центре нижней камеры повышается до 16% (по сравнению с 7% в сценариях а) и в)).

Рассмотрим энергетический спектр ионов на нижнем электроде. При отключенном ВЧ-напряжении данный электрод становится так называемым “плавающим” (“floating electrode”) – случай контакта проводящей незаземленной поверхности с плазмой. Около такой поверхности формируется электрическое поле, обеспечивающее равенство потоков электронов и ионов на поверхность. Соответствующую разность потенциалов между плавающим электродом и плазмой можно оценить, как [46]

(21)
${{V}_{{pl}}} - {{V}_{{fl}}} \cong \frac{1}{2}\frac{{k{{T}_{e}}}}{e}\ln \left( {\frac{{\pi {{M}_{i}}}}{{2{{m}_{e}}}}} \right).$

В случае плазмы аргона Mi = 40 а.е.м и ${{V}_{{fl}}} \cong 5.8k{{T}_{e}}{\text{/}}e$, где Te – температура электронов в области около приэлектродного слоя. Распределение потенциала около нижнего электрода показано на рис. 9. Из рисунка видно, что данное распределение соответствует случаю плавающего электрода с разностью потенциалов порядка 14 В (с учетом значения температуры электронов 2.4 эВ около электрода).

Рис. 9.

Распределение потенциала электрического поля на оси разряда около электрода.

Расчетный энергетический спектр ионов Ar+ представлен на рис. 10. Спектр имеет максимум на энергии, соответствующей разности потенциалов между плазмой и плавающим электродом. На рисунке виден вклад низкоэнергетичных ионов в спектр, так как, несмотря на маленькое давление газа 10 мТорр, ширина слоя около электрода составляет порядка 2.4 см. Ионы испытывают соударения во время движения к электроду.

Рис. 10.

Функция распределения ионов по энергии на плавающем электроде.

Для контроля энергетического спектра в экспериментальной системе на нижний электрод может подаваться также дополнительное переменное напряжение. Для описания такой системы фактически необходимо моделировать емкостной разряд в нижней камере. Такое моделирование требует хорошего разрешения по периоду приложенного ВЧ-напряжения. Стационарное решение может достигаться за тысячи периодов поля и процесс решения требует гораздо большего времени (часы и даже дни, в зависимости от размеров моделируемой системы). В дальнейшем настоящая модель индукционного разряда будет расширена на случай расчетов с ВЧ‑напряжением на нижнем электроде.

3.3. Результаты расчетов по температуре газа и зависимости от величины расхода газа

Многие модели индукционного разряда не включают в себя отдельных уравнений для температуры и протока газа: температура предполагается постоянной, а проток не учитывается, как, например, в [24, 28]. Тестовые расчеты в разд. 3.1 показали, что даже в таком приближении можно получить результаты, согласующиеся с экспериментальными данными. Однако включение в модель разряда уравнений для температуры и скорости газа позволяет аккуратно учесть влияние разогрева газа и времени пребывания газа в разряде на скорости химических реакций и параметры разряда. Это особенно важно при переходе от описания разрядов в инертном газе (в данной работе в аргоне) к описанию разрядов в смеси с молекулярными газами, такими как, например, кислород и хлор, активно используемыми в индукционных реакторах [47].

В данном разделе представлены результаты по температуре газа в зависимости от вкладываемой мощности и зависимость параметров плазмы от величины расхода газа. Все расчеты проводились в геометрии установки [11], рассмотренной в предыдущем разделе.

Существуют различные варианты граничного условия для температуры: фиксированная температура стенки, теплоизолированная стенка (равенство нулю теплового потока), заданная температура (например, комнатная) на определенном расстоянии от внешней границы камеры. Выбор конкретного варианта зависит от экспериментальных условий. В работе [11] температура стенки не была зафиксирована, стенка также не была теплоизолирована. В таких условиях для точного сравнения с экспериментом необходимы экспериментальные данные по измерению температуры, например, стенок камеры. Приблизительные измерения показывают, что нижняя камера нагревается до температуры порядка 400 К. Поэтому сравнение в табл. 3 приведено для распределения температуры в нижней камере, которое соответствует температуре 405 К на стенках (смотри рис. 7).

После сравнения с имеющимися экспериментальными данными, которое показывает применимость модели для описания данной экспериментальной установки, было также изучено изменение параметров разряда при увеличении вкладываемой в разряд мощности.

Зависимость температуры газа в центре разряда под катушкой от вложенной мощности в диапазоне от 100 до 1000 Вт представлена на рис. 11. Из рисунка видно, что при увеличении мощности возможен значительный разогрев газа: до температуры порядка 1000 К для вложенной мощности 1000 Вт. При условии заданного давления соответственным образом уменьшается плотность нейтрального газа, а при расчете в разряде аргона с молекулярными добавками возможно изменение роли химических процессов, зависящих от температуры.

Рис. 11.

Зависимость температуры газа в центре разряда от вложенной мощности.

Рис. 12.

Профиль концентрации плазмы на оси ($r = 0$) в зависимости от величины расхода газа Q = 0, 3, 20 sccm.

Также были проведены расчеты с различными значениями расхода газа: от значения 3 sccm, использованного в экспериментальных условиях [11] до 20 sccm.

При добавлении протока газа в модель давление в камере перестает быть постоянным в пространстве. Фактически, значение давления задается в месте, где газ выходит из камеры, а в остальном объеме рассчитывается исходя из системы уравнений. При этом появляется градиент давления в верхней камере. Для расхода 3 sccm увеличение давления в верхней части камеры составляет до 10% от базового уровня. Изменение профиля концентрации плазмы на оси разряда в зависимости от величины расхода газа Q = 0, 3, 20 sccm приведено на рис. 12. С ростом расхода плотность плазмы становится больше в верхней части камеры вследствие увеличения давления и плотности нейтрального газа. При этом в нижней камере плотность плазмы уменьшается, так как расчеты производятся при одной фиксированной мощности. Изменение других параметров плазмы: температуры электронов Te, потенциала плазмы Vpl, температуры газа Tgas, и плотности нейтрального газа N показано на рис. 13. Изменение данных параметров коррелирует с описанным выше градиентом давления. Учет протока газа особенно важен в случае моделирования разряда в смеси газов со сложной кинетикой нейтральных компонент (например, Ar/O2/Cl2), так как проток определяет время пребывания газа в разрядной области.

Рис. 13.

Профили: температуры электронов (а), потенциала плазмы Vpl (б), температуры газа Tgas (в) и плотности нейтрального газа N (г) на оси ($r = 0$) в зависимости от величины расхода газа Q = 0, 3, 20 sccm.

4. ЗАКЛЮЧЕНИЕ

В результате работы была построена и протестирована модель индукционного плазменного разряда в диффузионно-дрейфовом приближении. Проведены расчеты для разряда в аргоне для трех различных геометрий плазменных установок: простой цилиндрической геометрии [43], геометрии GEC [44] и геометрии с двумя камерами [11]. Результаты данных расчетов согласуются как с экспериментальными значениями, так и с другими численными моделями. В представленную в работе модель включены также уравнения для температуры и скорости нейтрального газа. Было показано, что учет этих физических процессов приближает результаты расчетов к экспериментально измеренным значениям в условиях сложной геометрии реактора [11].

Получены пространственные распределения плотности плазмы, температуры электронов и потенциала электрического поля. Проведен расчет температуры газа в зависимости от вложенной мощности. Приведено изменение параметров плазмы в зависимости от расхода газа.

Модель индукционного разряда дополнена кинетической моделью движения ионов в заданном электрическом поле с учетом столкновений с нейтральными частицами методом Монте-Карло. Показан энергетический спектр ионов Ar+ на плавающем электроде.

Построенная диффузионно-дрейфовая модель адекватно описывает плазменные характеристики и может быть использована при дальнейшем исследовании плазмы индукционного разряда для вычисления потоков интересующих компонент плазмы на подложку.

В дальнейшем планируется расширить модель на случай подачи ВЧ-напряжения на подложку, для контроля энергетического спектра ионов на обрабатываемой поверхности.

Исследование выполнено за счет гранта Российского научного фонда (РНФ № 18-72-00155).

Список литературы

  1. Hopwood J. // Plasma Sources Sci. Technol. 1992. V. 1 (2). P. 109.

  2. Laroussi M. // Plasma Process. Polym. 2005. V. 2 (5). P. 391.

  3. Lettry J., Aguglia D., Alessi J., Andersson1 P., Bertolo S., Briefi S., Butterworth A., Coutron Y., Dallocchio A., David N., Chaudetet E. et al. // Rev. Sci. Instrum. 2016. V. 87. P. 02B139.

  4. Fantz U., Heinemann B., Wünderlich D., Riedl R., Kraus W., Nocentini R., Bonomo F. // Rev. Sci. Instrum. 2016. V. 87. P. 02B307.

  5. Donnelly V.M., Kornblit A. // J. Vacuum Sci. Technol. A. 2013. V. 31. P. 050825.

  6. Athavale S.D., Economou D.J. // J. Vacuum Sci. Technol. B. 1996. V. 14. P. 3702.

  7. Shin H., Zhu W., Xu L., Donnelly V.M., Economou D.J. // Plasma Sources Sci. Technol. 2011. V. 20 P. 055001.

  8. Meichsner J., Wegner T. // European Phys. J. D. 2018. V. 72. P. 85.

  9. Adamovich I., Baalrud S.D., Bogaerts A., Bruggeman P.J., Cappelli M., Colombo V., Czarnetzki U., Ebert U., Eden J.G., Favia P., Graves D.B. et al. // J. Phys. D: A-ppl. Phys. 2017. V. 50. P. 323001.

  10. Александров А.Ф., Вавилин К.В., Кралькина Е.А., Павлов В.Б, Рухадзе А.А. // Физика плазмы. 2007. Т. 33. С. 816.

  11. Lopaev D.V., Rakhimova T.V., Rakhimov A.T., Zotovich A.I., Zyryanov S.M., Baklanov M.R. // J. Phys. D: Appl. Phys. 2018. V. 51. P. 02LT02.

  12. Uchida S., Takashima S., Hori M., Fukasawa M., Ohshima K., Nagahata K., Tatsumi T. // J. Appl. Phys. 2008. V. 103 (7). P. 073303.

  13. Baklanov M.R., de Marneffe J.-F., Shamiryan D., Urbanowicz A.M., Shi H., Rakhimova T.V., Huai Huang, Paul S. Ho // J. Appl. Phys. 2013. V. 113 (4). P. 041101.

  14. Nishida K., Mattei S., Mochizuki S., Lettry J., Hataya-ma A. // J. Appl. Phys. 2016. V. 119 P. 233302.

  15. Mattei S. Nishida K., Onai M., Lettry J., Tran M.Q., Hatayama A. // J. Computational Phys. 2017. V. 350 P. 891.

  16. Kawamura E., Birdsall C.K., Vahedi V. // Plasma Sources Sci. Technol. 2000. V. 9. P. 413.

  17. Makabe T. Advances in Low Temperature RF Plasmas: Basis for Process Design. Amsterdam: Elsevier, 2002.

  18. Birdsall C.K., Langdon A.B. Plasma Physics via Computer Simulation. N. Y.: McGraw-Hill, 1985.

  19. Tajima T. Computational Plasma Physics. Redwood city: Addison-Wesley, 1988.

  20. Cheng J., Ji L., Zhu Y., Shi Y. // J. Semicondactors. 2010. V. 31. P. 032004.

  21. Gogolides E., Sawin H.H. // J. Appl. Phys. 1992. V. 72. P. 3971.

  22. Hsu C.-C., Nierode M.A., Coburn J.W., Graves D.V. // J. Phys. D: Appl. Phys. 2006. V. 39. P. 3272.

  23. Kim H.C., Iza F., Yang S.S., Radmilović-Radjenović M., Lee J.K. // J. Phys. D: Appl. Phys. 2005. V. 38. P. R283.

  24. Brezmes A.O., Breitkopf C. // Vacuum. 2015. V. 116. P. 65.

  25. Lymberopoulos D.P., Economou D.J. // J. Res. National Institute Standards Technology. 1995. V. 100. P. 473.

  26. Райзер Ю.П. Физика газового разряда: учеб. руководство. М.: Наука. Гл. ред. физ.-мат. лит., 1987.

  27. Hagelaar G.J.M., Pitchford L.C. Solving the Boltzmann equation to obtain electron transport coefficients and rate for fluid coefficients models. Toulouse: Centre de Physique des Plasmas et de leurs Applications de Toulouse, 2005.

  28. Brezmes A.O., Breitkopf C. // Vacuum. 2014. V. 109. P. 52.

  29. Singh H., Graves D.B. // J. Appl. Phys. 2000. V. 88. P. 3889.

  30. Godyak V.A., Piejak R.B., Alexandrovich B.M. // Plasma Sources Sci. Technol. 2002. V. 11. P. 525.

  31. Mahoney L.J., Wendt A.E., Barrios E., Richards C.J., Shohet J.L. // J. Appl. Phys. 1994. V. 76. P. 2041.

  32. Lieberman M.A., Lichtenberg A.J. Principles of Plasma Discharges and Materia Processing. New Jercey: John Willey & Sons, Inc., 2005.

  33. Richards A.D., Thomson B.E., Sawin H.H. // Appl. Phys. Lett. 1987. V. 50. P. 492.

  34. Graves D.B., Kushner M.J., Gallagher J.W., Garscadden A., Oehrlein G.S., Phelps A.V. Database needs for modeling and simulation of plasma processing: National Research Council, Panel of Database Needs in Plasma Processing., 1996.

  35. Phelps A.V., Petrović Z.Lj. // Plasma Sources Sci. Technol. 1999. V. 8. P. R21.

  36. Dixon A.J., Harrison M.F.A., Smith A.C.H. // J. Phys. B. 1976. V. 9. P. 2617.

  37. Lee C., Lieberman M.A. // J. Vacuum Sci. Technol. A. 1995. V. 13. P. 368.

  38. Григорьян Г.М., Дятко Н.А., Кочетов И.В. // Физикаплазмы. 2015. Т. 41. С. 471.

  39. Фейнман Р., Лейтон Р., Сэндс М. Фейнмановские лекции по физике. Вып. 6. Электродинамика. М.: Мир, 1964.

  40. www.grc.nasa.gov/www/CEAWeb/.

  41. COMSOL 5.2, Userbook.

  42. Phelps A.V. // J. Appl. Phys. 1994. V. 76 P. 747.

  43. Godyak V.A. // Electron Kinetics and Applications of Glow Discharges /Eds. U. Kortshagen, L. D. Tsendin. N. Y.: Plenum Press, 1998. P. 241.

  44. Miller P.A., Hebner G.A., Greenberg K.E., Pochan P.D., Aragon B.P. // J. Res. National Institute Standards Technology. 1995. V. 100 P. 427.

  45. Zotovich A.I., Rezvanov A., Chanson R., Zhang L., Hacker N., Kurchikov K.A., Klimin S., Zyryanov S.M., Lopaev D.V., Gornev E., Clemente I., Miakonkikh A., Maslakov K.I. // J. Phys. D: Appl. Phys. 2018. V. 51 P. 325202.

  46. Boswell R.W., Morey I.J. // Appl. Phys. Lett. 1988. V. 52 P. 21.

  47. Bogart K.H.A., Donnelly V.M. // J. Appl. Phys. 1999. V. 86. P. 1822.

Дополнительные материалы отсутствуют.