Адаптивное изменение массы модельных частиц при моделировании тлеющего ВЧ-разряда в силановой плазме Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 13, № 1, 2008

Адаптивное изменение массы модельных частиц при моделировании тлеющего ВЧ-разряда

*.9 _ *

в силановои плазме*

В. А. Вшивков, Г. Г. Лазарева, А. В. Снытников Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: lazareva@ssd.sscc.ru

A method is proposed that enables to simulate the recombination of positive and negative ions in an accurate way, even at low densities. The recombination in a cell is done by changing the mass of superparticles in this cell. Thus it is possible to recombine any desired amount of real particles even if it is less than one superparticle. Moreover, the number of superparticles in each cell is kept nearly constant by adding zero-mass superparticles (if a cell has too few superparticles) or by removing a superparticle with the increase of mass for each in the rest of superparticles (if a cell has too many superparticles). It is shown that using this method a simulation of the ion recombination results in a lower level of non-physical fluctuations. The results of computations are compared to the experimental data and the results of hydrodynamical simulations.

Введение

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

Численные модели, используемые для исследования низкотемпературной (в частности, силановой) плазмы могут быть разделены на три группы по методу решения кинетического уравнения Больцмана: гидродинамические, кинетические и гибридные. При кинетическом подходе [1,2] вычисляются зависящие от времени и координат решения уравнения Больцмана, которые дают распределения скорости электронов и ионов. Уравнение Больцмана решается либо непосредственно (с помощью так называемого двучленного приближения [3]), либо с использованием статистических подходов (метод Монте-Карло). При том, что кинетические модели требуют больших вычислительных затрат, они в меньшей степени зависят от изначальных предположений и дают более точные результаты.

*Работа выполнена при частичной финансовой поддержке гранта NWO-Plasma (Dutch-Russian NWO-GRID project, contract NWO-RFBS 047.016.007, Dutch-Russian NWO-Plasma project, contract NWO-RFBS 047.016.018), программы СО РАН по суперЭВМ, программы Рособразования “Развитие научного потенциала ВШ” (проект РНП.2.2.1.1.3653, проект РФФИ 05-01-00665).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

1. Постановка задачи

Рассмотрим кинетическое уравнение Больцмана, дополненное уравнением Пуассона для потенциала в цилиндрической системе координат:

f + vV/„ + ^ f = St{/„},

dt m dv

1

r dr \ dr ) + r2

1 д 2Ф 5 2Ф

+ = 4пр,

(1)

dz2

где р — плотность, v — скорость, E = —VФ — напряженность электрического поля, Ф — потенциал, /а — функция распределения частиц сорта a, St{f} — столкновительный член.

На рис. 1 изображена расчетная область. Сверху — ВЧ-электрод, на который подается потенциал вида

U(t) = Urf cos ut,

в настоящей работе 20 B < URF < 70 В, и = 58 МГц (частота ВЧ-поля). Снизу — заземленный электрод. На стенках задано условие равенства нулю нормальной производной потенциала. Электроны и ионы могут покидать область через электроды, от боковых стенок модельные частицы отражаются, на электродах — поглощаются. Размер области по радиусу Rm , по z-координате — Zm. Для расчета были взяты следующие параметры: Rm = 10 см, Zm = 2.5 см. В этой области введена цилиндрическая сетка Nr х Nv х NZ = 50 х 32 х 30. Расчеты проводились на кластере НКС-160 (Сибирский суперкомпьютерный центр в Новосибирске).

Рис. 1. Конфигурация расчетной области

2. Модель силановой плазмы ВЧ-разряда

В работе используется физическая модель тлеющего ВЧ-разряда в силановой плазме, представленная в [4, 5]. Более подробное описание физико-химических процессов в силановой плазме можно найти, например, в [6] Рассматриваемая модель силановой плазмы включает в себя ионы ВіН+, БіИ-, Н+ и электроны. Их движение происходит на фоне нейтрального газа, представляющего собой смесь силана БіН4 и водорода Н2. Концентрация нейтральных частиц (молекул силана и водорода) считается постоянной, их движение в рамках данной работы не рассматривается. Плотность заряженных

частиц на шесть порядков меньше плотности нейтрального газа, поэтому рассматриваются только столкновения электронов с молекулами газа. Давление газовой смеси 16 Па (0.00016 атм), суммарная концентрация силана и водорода 1016 см-3, концентрация электронов 108 см-3, суммарная концентрация ионов силана и водорода 1010 см-3. Оценка длины свободного пробега в нейтральном газе для электронов 0.1 см, для ионов

0.0001 см.

В модели рассматриваются следующие процессы:

— упругое столкновение e- + SiH4 ^ SiH4 + e-;

— ионизация e- + SiH4 ^ SiH+ + 2H + 2e-;

— колебательное возбуждение e- + SiH4 ^ SiHg + e-;

— прилипание e- + SiH4 ^ SiH- + H;

— диссоциация e- + SiH4 ^ SiH+ + 2H + 2e-;

— упругое столкновение e- + H2 ^ H2 + e-;

— ионизация e- + H2 ^ H+ + 2e-;

— колебательное возбуждение e- + H2 ^ Hg + e-;

— диссоциация e- + H2 ^ 2H + e-.

Для силана существует две реакции колебательного возбуждения, различающиеся порогом реакции, а для водорода — три. Образовавшиеся в реакциях радикалы (атомы водорода) в реальном разряде поглощаются на стенках камеры. Так как радикалы не имеют заряда, они не учитываются при расчете плотности заряда и потенциала. Таким образом, новые модельные частицы образуются только в результате реакции ионизации. Во всех реакциях, кроме упругого столкновения, электрон теряет энергию, равную порогу реакции. В настоящей работе модельными частицами представлены следующие сорта реальных частиц: электроны и ионы SiH+, SiH-, H+. В каждой ячейке изначально находится 100 модельных частиц каждого сорта, всего 19.2 млн частиц.

Моделирование рассеяния электронов на атомах водорода и силана во всех реакциях производится с помощью метода нулевых столкновений (null collision technique [1, 2]). При этом используются известные из эксперимента сечения [6] столкновения и пороги реакций. Кроме того, модель включает в себя рекомбинацию положительных и отрицательных ионов:

SiH+ + SiH- ^ SiH2 + SiH3,

H+ + SiH- ^ H2 + SiH3. (2)

Заряженные частицы (положительные ионы и электроны) поглощаются на стенках камеры. Отрицательные ионы не могут приблизиться к электродам в силу наличия приэлектродного слоя. Поэтому для моделирования стационарного разряда необходимо учитывать рекомбинацию как единственную возможность сокращения числа отрицательных ионов.

Расчет должен продолжаться до тех пор, пока разряд не станет стационарным (пока средний по времени профиль концентрации ионов и электронов не перестанет меняться). Так как масса иона намного (для иона в 55 080 раз) больше массы электрона, то при задании времени расчета нужно ориентироваться именно на ионный компонент плазмы. Таким образом, время расчета должно значительно (в несколько десятков раз) превышать время пролета иона через межэлектродный промежуток, при указанных выше параметрах расчета это более 1000 периодов ВЧ-поля, или 2 х 10-5 с.

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

более одной ячейки за временной шаг (это необходимо для устойчивости расчета с применением метода частиц-в-ячейках). В результате число временных шагов в расчете оказывается очень большим (более миллиона). Эта трудность преодолевается в настоящей работе проведением расчетов на суперЭВМ.

3. Описание метода

Приведем описание новой модификации метода частиц, основанной на адаптивном изменении масс модельных частиц.

Распределение частиц каждого сорта (при Ь = 0) выполняется следующим образом. Действительных частиц в каждой ячейке N штук. Число N выбирается исходя из оценки уровня численных шумов N-1/2. В настоящей работе N = 100, что обеспечивает уровень численных шумов ~ 10 %. Масса и заряд пропорциональны плотности данного вещества в ячейке. Задается гауссово распределение частиц по скоростям с температурой, изначально равной температуре нейтрального газа (от 500 до 1000 °С).

Отношение р = — - константа для данного сорта частиц. Частица хранит или массу т, т

или заряд q.

Если частиц данного сорта в ячейке нет, но их появление возможно (например, за счет ионизации), то распределяются фиктивные частицы также в количестве N. Масса и заряд фиктивных частиц равны нулю. Распределение фиктивных частиц по скоростям должно соответствовать распределению частиц, из которых возникают данные фиктивные частицы, т. е. средняя скорость должна оставаться неизменной при добавке новых частиц. Движение фиктивных частиц происходит по инерции, с зеркальным отражением от стенок.

Фиктивные частицы нужны для того, чтобы поддерживать постоянным (около 100) число частиц в ячейке, а именно в тех ячейках, в которых реальных частиц соответствующего сорта на данный момент нет. Если в такую ячейку попадает реальная частица (одна), то плотность в этой ячейке будет вычисляться с большими погрешностями (величина численных шумов около 100 %). При наличии фиктивных частиц в этой ячейке они превращаются в реальные, и вычисление плотности происходит с более низким уровнем численных шумов (10 %).

Если же каждый раз, как только частиц в данной ячейке оказывается мало, просто добавлять частицы (доводить их число до нормы — до 100), то число частиц в программе (требуемый объем памяти) будет неограниченно расти.

Если в ячейке слишком мало действительных или фиктивных частиц

Nt <

то необходимо добавлять частицы. Здесь N — текущее количество частиц в ячейке,

— модельный коэффициент.

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

Приведем пошаговое описание алгоритма.

1. Вычисляется полная масса в ячейке.

2. Вычисляется средняя скорость.

3. Добавляются новые частицы ^ — N (Ь) штук) с необходимым разбросом по скоростям (средняя скорость равна нулю).

4. Вычисляется масса одной частицы в ячейке делением полной массы на полное количество частиц (и старых, и новых). Эта масса присваивается каждой частице в ячейке.

5. К случайным скоростям новых частиц добавляется средняя скорость частиц в ячейке.

Фиктивные частицы добавляются так же, как в начальный момент.

Если, наоборот, в ячейке слишком много действительных или фиктивных частиц

^ > кЛе\ .М,

то необходимо удалять частицы. Здесь Nt — текущее количество частиц в ячейке, к^е\ — модельный коэффициент. Удаление частиц происходит аналогично добавке с сохранением массы (заряда) и средней скорости частиц в ячейке, только вместо добавки в пункте 3 алгоритма удаляется kde\N — N частиц. Значения коэффициентов подбираются так, чтобы и удаление, и добавка частиц происходили как можно реже, например, при к^е\ = 2, кааа = 0.5 частицы добавляются в среднем в 9 % ячеек, а удаляются в 3 % ячеек в течение одного временного шага.

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

4. Моделирование процесса рекомбинации

В разделе 3 дается принципиальное описание алгоритма адаптивного изменения массы модельных частиц. Далее приведем описание его работы в конкретной ситуации (рекомбинация ионов).

Рекомбинация положительных ионов с концентрацией п+ и отрицательных ионов с концентрацией п- описывается уравнением

пп -

~д~ = кГесП+п-. (3)

Здесь кгес — коэффициент рекомбинации (для реакции + БШ- коэффициент ре-

комбинации кгес = 1.2 • 10-7 см3/с). Из этого уравнения можно определить количество ионов, которые исчезают в результате рекомбинации в течение небольшого промежутка времени ДЬ в малом объеме ДУ:

Дп+ = кгесп+п- ДЬДУ. (4)

В рамках рассматриваемой модели возможна рекомбинация отрицательных ионов

БШ- с положительными ионами БШ+Г и Н+. Рассмотрим рекомбинацию ионов БШ- и

81Н+.

Трудность моделирования рекомбинации заключается в переносе этого процесса на модельные частицы. Это можно проиллюстрировать следующим образом. Если в формуле (4) интервал времени соответствует временноому шагу процесса моделирования ДЬ = т = 1.68 • 10-10 с, а малый объем — объему ячейки сетки ДУ = ткгк1ркг 4• 10-3 см3, то при характерной плотности ионов п+ = п- = 1010 см-3 число рекомбинирующих ионов в одной ячейке оказывается равным примерно 10. Это, очевидно, намного меньше

числа реальных ионов, соответствующих одной модельной частице (при указанных в разделе 1 параметрах расчета это примерно 106 ионов). Таким образом, если в ячейке содержится 100 тыс. модельных частиц (что нереально), то рекомбинация моделируется удалением одной из них на каждом временном шаге (если скорость рекомбинации постоянна). Однако если частиц меньшее количество, то необходимо распространять рекомбинацию или на большее количество ячеек, или на большее количество временных шагов. При указанных параметрах расчета (48 тыс. ячеек сетки, 10 рекомбинирующих реальных ионов в ячейке, 106 реальных ионов на одну модельную частицу) для моделирования рекомбинации должна удаляться одна модельная частица во всей области в течение двух временных шагов.

Таким образом, в одной из ячеек сетки будет удалена частица, что приведет к значительному изменению плотности в этой ячейке, в то время как в окружающих ячейках плотность не изменится. Это означает, что будут внесены дополнительные нефизические флуктуации плотности (с точки зрения физики темпы изменения плотности в ре-

Ан

зультате рекомбинации, характеризующиеся величиной-----в соседних ячейках должны

н

быть близкими по величине).

Более того, изменение плотности в той ячейке, в которой была удалена частица, будет значительно (в приведенном примере, в 105 раз) больше, чем должно быть в результате именно рекомбинации.

Описанный в разделе 3 подход позволяет моделировать рекомбинацию ионов более естественным образом: появляется возможность удалить любое малое количество реальных частиц, даже если это меньше, чем соответствует одной модельной частице. Это происходит следующим образом.

Плотность ионов в ячейке уменьшается на величину, равную Ан+. Это означает, что масса каждой из модельных частиц, принадлежащих данной ячейке, умножается

Рис. 2. Максимальное изменение плотности положительных ионов, нормированное на величину плотности в зависимости от номера временного шага, х(Ьп): линия с маркерами — классический метод частиц-в-ячейках, все модельные частицы имеют равные фиксированные массы; гладкая линия — алгоритм с адаптивным изменением масс модельных частиц. В первом случае рекомбинация проводится только на каждом десятом временном шаге

на коэффициент

к

1 — Дп+

п

+

В этом случае отношение изменения плотности к величине плотности будет одинаковым по всем ячейкам с номерами г,р,к:

(Дп+)

г,р,к

(п+)г,р,к (п )г,р,к

кГе.С.ДЪ'

(5)

На рис. 2 показаны максимальное значение величины, характеризующей нефизическое изменение плотности в результате рекомбинации:

х(гп)

тах

г,р,к

(Дп+)

)г,р,к

п

+

г,р,к

г = 1,..., р = 1

N, к

1,

N

z.

(6)

Здесь Ьп — номер временного шага. Изменение плотности в выражении (6) является результатом только рекомбинации, другие процессы, в частности, ионизация и перенос вещества, не учитываются. Таким образом, если удаление одной модельной частицы в некоторой ячейке сетки соответствует рекомбинации в объеме нескольких ячеек в течение нескольких временных шагов, то величина £ будет существенно больше реального физического значения изменения плотности, определяемого из выражения (3) (рис. 2, линия 1). Наоборот, в том случае, когда рекомбинация проводится отдельно в каждой ячейке сетки, £ примерно равно физическому изменению плотности кгесДЬп+ ~ 10-7, как видно из рис. 2, линия 2.

5. Результаты

На рис. 3 дана функция распределения электронов по энергиям. Функция распределения f (е) показывает число электронов, имеющих энергию в интервале [е, е + de], энергия задается в электронвольтах. В настоящей работе функция распределения используется только для сравнения результатов моделирования с другими работами, в частности, с гидродинамическим моделированием. Сравнение производится по отношению максимального и минимального значений и по ширине максимума.

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

I, нА 600

400

200

10

20

А 1 • 2

Р, Вт

Рис. 4. Зависимость тока ионов на заземленный электрод от мощности для частоты ВЧ-поля 58 МГц в расчетах (треугольники) в сравнении с экспериментом [7] (круги)

На рис. 3 видно, что максимальное количество электронов находится в интервале 1.5-4 эВ, энергия электронов в максимуме функции распределения электронов по энергиям равна 2 эВ. Распределение электронов по энергиям приближенно соответствует результатам гидродинамического моделирования, представленным в работе [5], в которой функция распределения электронов по энергиям вычисляется путем приближенного решения уравнения Больцмана с помощью двучленного приближения.

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

Результаты этих расчетов приведены на рис. 4 в сравнении с результатами эксперимента [7]. На рис. 4 видно, что характер зависимости тока от мощности в расчетах тот же, что и в эксперименте. Различия в данных отчасти объясняются неполным соответствием геометрии расчетной области и геометрии экспериментальной установки.

Вывод

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

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

[1] BlRDSALL C.K. Particle-in-Cell Charged-Particle Simulations Plus Monte Carlo Collisions With Neutral Atoms, PIC-MCC // IEEE Trans. Plasma Sci. 1991. Vol. 19, N 2. P. 65-83.

[2] Vahedi V., SuRENDRA M. A Monte Carlo collision model for the particle-in-cell method: applications to argon and oxygen discharges // Computer Physics Communications. 1995. Vol. 87. P. 179-198.

[3] ШвЕйГЕРТ В.А. Численное моделирование стационарной функции распределения электронов в слабоионизованном газе в неоднородных электрических полях // ПМТФ. 1989. № 5. C. 3-7.

[4] Nienhuis G.J., Goedheer W.J., Hamers E.A.G., Sark W.G.J.H.M. van, Bezemer J. A self-consistent fluid model for radio-frequency discharges in SiH4-H2 compared to experiment // J. Appl. Phys. 1997. Vol. 82. P. 2060.

[5] Nienhuis G.J. Plasma models for silicon deposition. Ph.D. Thesis. Utrecht Univ., 1998.

[6] Perrin J., Leroy C., Bordage M.C. Contrib // Plasma Phys. 1996. Vol. 36. P. 3.

[7] Абрамов А.С., Виноградов А.Я., Косарев А.И., Смирнов А.С., Орлов К.Е., Шутов М.В. Исследование ионной бомбардировки пленок аморфного кремния в процессе

плазмохимического осаждения в высокочастотном разряде // ЖТФ. 1998. Т. 68, № 2.

С. 52-59.

Поступила в редакцию 23 марта 2007 г., в переработанном виде —13 июня 2007 г.