CC BY
112
В. М. Грабов, Е. В. Демидов, Е. В. Усынин
МОДЕЛИРОВАНИЕ ПРОЦЕССА РОСТА ПЛЕНОК ВИСМУТА НА ПОДЛОЖКЕ ИЗ СЛЮДЫ
Работа выполнена при финансовой поддержке Министерства образования и науки РФ в рамках реализации Аналитической ведомственной целевой программы «Развитие научного потенциала
высшей школы (2009-2010 годы)» (грант № 2.1.1/3847) и Федеральной целевой программы «Научные и научно-педагогические кадры инновационной России»
(гос. контракт от 22 марта 2010 г. № 02.740.11.0544).
Приведено описание модели, реализованной с использованием метода Монте-Карло, позволяющей выявить закономерности влияния параметров технологического режима получения пленки висмута на форму, размеры и кристаллографическую ориентацию кристаллитов. Описаны особенности реализации и оптимизации модели. Приведено сравнение полученных результатов моделирования с особенностями роста реальных пленок висмута на слюде, получаемых методом термического испарения.
Ключевые слова: моделирование, висмут, слюда, эпитаксия, метод Монте-Карло, тонкие пленки, кристаллиты, адсорбция, диффузия.
V. Grabov, E. Demidov, E. Usynin
SIMULATION OF BISMUTH THIN FILMS GROWTH PROCESS ON MICA
The paper describes a Monte-Carlo model that allows to reveal natural laws of the impact of the technological conditions parameters of the production of bismuth thin films on the bismuth crystallite shape, size and crystallographic orientation. The model implementation and optimization features are also described. A comparison is given of the model simulation results and the natural bismuth thin films produced by the method of thermal spraying on mica.
Key words: simulation, bismuth, mica, epitaxy, Monte-Carlo method, thin films, crystallo-graphic orientation, adsorption, diffusion.
Введение
Развитие электроники способствует непрерывному прогрессу и углублению научных знаний в области пленочного материаловедения. В то же время усложнение структуры функциональных устройств микроэлектроники на основе эпитаксиальных пленок предъявляет жесткие требования к методам получения, толщинам и свойствам отдельных слоев. Оптимизации условий выращивания тонких пленок способствуют все более широко применяемые методы математического моделирования процессов роста пленок, использующие возможности современных быстродействующих ЭВМ.
В физике конденсированного состояния широко применяются два метода моделирования: Монте-Карло и молекулярной динамики. Нерешеточные модели успешно изучаются обоими методами, процесс роста кристаллов — в основном решеточным вариантом метода Монте-Карло [1].
Использование аппарата теории вероятности для решения прикладных задач с помощью ЭВМ приобрело большую популярность. Метод Монте-Карло используют для изучения различных аспектов процесса роста тонких пленок. Изучение роли одномерной диффузии приводится в работе [2], комплексный анализ использования метода Монте-Карло
при моделировании роста полупроводниковых пленок можно найти в работе [1]. В данной работе представлена реализованная методом Монте-Карло модель процессов, происходящих во время выращивания пленки висмута методом термического напыления.
Характерной особенностью пленок висмута, получаемых методом термического испарения в вакууме на кристаллических подложках из слюды, является ориентирующее действие слюды на пленку висмута: ось С3 в кристаллитах висмута ориентирована перпендикулярно подложке, при этом существуют 2 энергетически наиболее выгодные кристаллографические ориентации кристаллита в плоскости подложки. Эти 2 ориентации отличаются противоположным направлением кристаллографической оси С2 [3].
Существующие программные комплексы для моделирования роста пленок не учитывают указанные особенности процесса формирования пленок висмута на слюде, в связи с чем не отражают в полной мере структуру пленки висмута.
Исследование структуры пленок показало, что размеры кристаллитов, из которых состоит пленка, могут варьироваться в большом интервале в зависимости от температуры подложки, толщины пленки, температуры отжига и изменяются от величины, примерно равной толщине пленки, до размеров, на порядок превышающих ее толщину [4]. Кристаллиты пленок, полученных при оптимальных условиях [4], представляют собой переплетение блоков двух ориентаций, имеющих большую длину и сложную конфигурацию. Появляется необходимость создания модели, описывающей рост таких пленок. Главной задачей модели является описание влияния параметров технологического режима получения пленки висмута на форму, размеры и кристаллографическую ориентацию кристаллитов.
Разработанный программный комплекс позволяет анализировать процесс роста пленок висмута в послойном режиме, находить оптимальные условия получения тонких пленок висмута без использования физического эксперимента, что, в свою очередь, значительно ускорит и удешевит проведение экспериментов.
1. Описание модели
Как известно, процесс роста тонких пленок является фазовым переходом первого рода. Кинетика фазовых переходов первого рода — сложный многостадийный процесс, кото -рый сопровождают различные нелинейные явления. К таким стадиям обычно относят стадии зародышеобразования, сепаратного роста зародышей новой фазы, коалесценцию и позднюю стадию, т. е. оствальдовское созревание, на которой рост более крупных островков осуществляется за счет растворения более мелких. Указанные процессы имеют различные масштабы времени. Наиболее быстро протекает стадия зарождения, затем стадия сепаратного роста и т. д. Эта иерархия времен означает, что быстрые процессы успевают «подстроиться» под медленные [9]. По существу, моделирование быстрых процессов дает начальные условия для моделирования более медленных процессов. Как быстрые, так и медленные процессы при использовании аналитического моделирования описываются совокупностью интегро-дифференциальных уравнений.
Применение метода Монте-Карло позволяет решать задачи кинетики фазовых переходов с большим числом переменных в неоднородной многочастичной системе. Существенно, что применение метода Монте-Карло не требует формулировки интегро-диф-ференциальных уравнений, описывающих многочастичную систему. На практике формулировка подобных уравнений часто представляет значительные трудности, а с помощью метода Монте-Карло можно решать задачи без аналитического исследования вероятностной модели [7].
При построении с помощью метода Монте-Карло модели роста тонких пленок и преобразования рельефа на атомном уровне формулируются типы событий, которые могут происходить. К таким событиям относятся адсорбция частиц на поверхность, диффузионный скачок частицы на свободное место, реиспарение частицы и т. д. Для корректной работы модели необходимо знать вероятность того или иного события.
В соответствии с принципом Больцмана, любая термодинамическая система занимает
возможное состояние с энергией Е с вероятностью, пропорциональной ехр(_Ег- / ^ьТ). Выражение Гиббса для вероятности Рг- занять системой состояние с энергией Е имеет следующий вид:
Р ОехР(~Е /кьТ) (1.
1 2 ’
где О — число различных состояний системы с одинаковой энергией Е, кь — постоянная Больцмана, Т — температура подложки, 2 — статистическая сумма.
Из вышесказанного следует, что для построения модели необходимо знать энергии взаимодействия структурных элементов в различных состояниях. В их числе поверхностная энергия границы раздела висмут — слюда, поверхностная энергия границы раздела кристаллитов пленки и модельное изменение этих энергий в зависимости от угла поворота кристаллографической оси С2 относительно направления преимущественной ориентации, задаваемого подложкой.
1.1. Расчет энергии взаимодействия поверхностей висмут — слюда
Трехмерный островок с контактным углом р находится в равновесии с подложкой и газообразной фазой при условии равенства сил, действующих на точку соприкосновения твердой, адсорбированной и газообразной фаз (соотношения Юнга):
у* = ус С0Бр + ус^, (2)
где у8 — поверхностная энергия подложки, ус — поверхностная энергия осаждаемого материала, ус _* — поверхностная энергия границы раздела адсорбат-подложка, р — угол
наклона грани островка к подложке [6].
В расчетах использовались следующие значения поверхностной энергии: для висмута уй = 3,88 х 10-1 Дж/м2, для слюды у* = 7,8х10-1 Дж/м2 [11]. Угол р вычислялся на основе
экспериментальных данных, полученных с использованием метода атомно-силовой микроскопии (АСМ).
На рис. 1а представлено АСМ-изображение островковой пленки висмута. При этом почти вся подложка покрыта сросшимися островками висмута, однако присутствуют участки, не содержащие островков. Исходя из АСМ-изображения границы подложка-островок (рис. 1б), был определен контактный угол островка р & 14°, что в соответствии с выражением (2) дает значение ус_* = 3,5 х 10-1 Дж/м2. С учетом количества атомов в слое значения энергии ус,у*,ус_* можно выразить в эВ/атом.
Рис. 1а. АСМ-изображение участка островковой пленки висмута
Рис. 1б. Изменение высоты среза пленки по линии, выделенной на рис. 1а
1.2. Описание изменения энергии взаимодействия висмут — слюда в зависимости от кристаллографической ориентации зародыша относительно подложки
При возникновении на подложке зародыша критического размера он может иметь различную кристаллографическую ориентацию относительно подложки (угол поворота оси С2 относительно энергетически наиболее выгодного направления). Каждая кристаллографическая ориентация характеризуется своей поверхностной энергией границы раздела. Из всего многообразия кристаллографических ориентаций осей С2 вследствие тригональ-ной симметрии наиболее вероятными являются лишь 6 из них. Зародыши с одной из этих 6 ориентаций имеют наименьшую энергию границы раздела зародыш-подложка, следовательно, и наибольшую вероятность возникновения. Вследствие тригональной симметрии ориентации, трансформированные поворотом вокруг оси С3 на 120o, кристаллографически эквивалентны. Поэтому в модели учитываются 2 наиболее вероятные ориентации. Они выделены пунктиром (рис. 2). В дальнейшем одна из 2 наиболее вероятных ориентаций оси С2 (не важно какая именно) берется в качестве начала отсчета, т. е. термин «угол поворота оси С2» какого-либо островка или зародыша будет обозначать угол между этой энергетически наиболее выгодной ориентацией оси С2 и другой произвольной ориентацией оси С2 островка или зародыша.
Поверхностную энергию границы раздела зародыш — подложка Es в зависимости от кристаллографической ориентации зародыша относительно подложки (угла поворота оси С2) можно описать уравнением
ES = Esa + Er (T,a) (3)
где Esa =Yd_s = 3,5 x 10 1 Дж/м2 — изотропная составляющая поверхностной энергии границы раздела зародыш — подложка,
Er (T,a) = (k • T + b) • f (a), (4)
— анизотропная составляющая поверхностной энергии границы раздела зародыш — подложка, зависящая от температуры подложки T и угла a поворота оси С2. Коэффициенты k, b
линейного множителя в выражении (4) подбираются, исходя из экспериментальных данных.
Рис. 2. Наиболее вероятные кристаллографические ориентации оси С2 фигур роста
Функция/ (а), описывающая модель изменения поверхностной энергии границы раздела зародыш — подложка в зависимости от угла поворота кристаллографической оси С2 зародыша относительно задаваемого подложкой, должна удовлетворять следующим условиям:
- функция/(а) должна быть определена на всем интервале а е [0; 2п];
- / (а) должна иметь локальные минимумы в точках а = 0 и а = п;
- область значений функции /(а) должна лежать в интервале [0, 1].
В процессе моделирования использовались следующие функции, удовлетворяющие представленным выше условиям:
2 п
—а,0 <а< — п 2
2 „ п
---а + 2,— <а<п
= < п 2
' 2 „ 3п
—а- 2,п<а<— п 2
2 3п
---а + 4,— <а < 2п;
_ п 2
= Нп(а)|;
[0,а = 0 V а = п
[1,а Ф 0 ла Ф п.
Расчет вероятности возникновения на поверхности подложки зародыша с углом поворота оси С2 (а) осуществляется с помощью выражения Гиббса (1), в котором используется энергия границы раздела зародыш — подложка (Е).
График плотности вероятности возникновения на поверхности подложки зародыша с углом поворота оси С2 имеет следующий вид (рис. 3). Пики на графике изображают наиболее вероятные углы поворота оси С2. Высота пиков определяется коэффициентами к и Ь
/ (а) / (а) / (а)
эВ 1 эВ
в (4). На графике представлена функция при к = -0,00265(-------------К ), Ь = 1,44(------), тем-
атом
атом
пературе подложки 100°С и /(а) = |з£п(а)|.
- 1
0,025 ^ Г
0,020 -
Г"
1
1
0 015 1
щ
0 010 1
1
1
: 1__
' ^“
1 1 1 1
Рис. 3. График плотности вероятности угла поворота оси С2
1.3. Описание изменения поверхностной энергии границы раздела островков висмута в зависимости от угла поворота кристаллографической оси С2
Процесс роста пленки висмута при взаимодействии со слюдой происходит только на самом раннем этапе роста пленки, до образования сплошного слоя. После образования сплошного слоя в модели используется только энергия взаимодействия атомов висмута с висмутом.
В процессе роста островков висмута они начинают соприкасаться друг с другом, в результате чего у них появляются общие границы раздела. Эта граница характеризуется поверхностной энергией, зависящей от углов поворота кристаллографической оси С2 каждого из островков. Поверхностную энергию границы раздела одного островка можно представить уравнением:
Е,
'ЕУ
3=1
N
зт(-
(5)
2
где Ееу — энергия связи в кристалле висмута (2,15 эВ/атом) [8], Ытах — максимально возможное количество соседних островков, а — угол поворота оси С2 островка, для которого рассчитывается вероятность, а}- — угол поворота оси С2 каждого конкретного соседнего
островка, суммирование ведется по всем соседним островкам.
В процессе роста островков, когда подложка еще не покрылась сплошным слоем висмута, а границы раздела между островками уже начали появляться, расчет итоговой поверхностной энергии производится суммированием поверхностной энергии границы раздела подложка — адсорбат, рассчитанной по формуле (3), и поверхностной энергии границы раздела между островками висмута (5). Таким образом, итоговая поверхностная энергия в данном случае равна:
Е = Е + Ем.
40
1.4. Описание диффузии
Модель диффузии состоит из нескольких этапов. Первый этап состоит в определении вероятности диффузионного скачка. Для этого необходимо вычислить энергию связи диффундирующего атома со всеми соседними атомами:
Е° = х ы->™' ’ <7)
тах
где Ео — энергия, необходимая для совершения диффузионного скачка, ЕЕу — энергия связи атомов висмута <2,15 эВ/атом) [8], Лгтах — максимально возможное количество соседних атомов, ЫпеЛоиг — реальное количество соседних атомов. Далее вероятность диффузии находим по формуле
PD _ Kosc exp(-Ed /kbT\ (8)
где Kose — связанный с особенностью реализации модели коэффициент, равный количеству тепловых колебаний атомов за время между двумя последовательными актами адсорб-
UOSC ^ .. Л „„12 — 1ч.-,--.
ции. ^osc - h х W х V , где uOsC — частота тепловых колебаний решетки (^ 10 сек ) [6],
H — ширина модельной решетки в атомоместах, W — длина модельной решетки в атомо-местах, V — скорость осаждения висмута (монослоев/сек).
Следующий шаг состоит в нахождении направления диффузионного скачка. Направление скачка определяется по принципу минимизации свободной поверхностной энергии. Для каждого возможного нового положения атома рассчитывается его свободная поверхностная энергия в этом положении, исходя из количества соседних атомов и количества сво-
eev
бодных связей, Ef _ N х Nvacant, где Ер — свободная поверхностная энергия, EEv —
max
энергия связи атомов висмута, Nmax — максимально возможное количество соседних атомов, Nvacant — число свободных связей. Далее для определения вероятности того или иного направления диффузии используется все та же формула (1).
Последний шаг заключается в определении принадлежности диффундировавшего атома к зародышу той или иной кристаллографической ориентации.
2. Реализация модели
Для моделирования стохастического процесса роста тонких пленок применялся метод Монте-Карло. Первоначально задаваемыми параметрами модели являются размер моделируемой пленки, температура подложки, скорость роста и энергии связи атомов висмута с кристаллами висмута и слюды. Создается модель абсолютно гладкой подложки размером N х N. Для имитации бесконечной плоскости задаются периодические граничные условия. Для зародыша с координатой <1, 1) соседним зародышем слева будет <1, ^, а сверху — (№, 1). Выбор места события на поверхности и определение вида события проводились с помощью датчика псевдослучайных чисел, генерирующего случайные числа в интервале 0-1. Два случайных числа выбираются для определения координат места события. Для по-
лучения координаты ячейки, в которой разыгрывается одно из элементарных событий, каждое случайное число умножается на число мест по осям.
Как следует из эксперимента [5], при термическом напылении пленок висмута на подложки из слюды формируются блоки двух устойчивых ориентаций с противоположным направлением осей С2. В модели островки, формирующие эти блоки, обозначаются черным и белым цветом. Кроме того, существует некоторое количество островков, формирующих кристаллиты с неустойчивой ориентацией кристаллографической оси С2, они обозначаются в модели серым цветом.
Следующий шаг при моделировании — адсорбция атомов. В процессе роста пленок висмута методом термического напыления в вакууме атомы адсорбируются из пересыщенного пара висмута [10]. Поэтому акты реиспарения атомов практически отсутствуют. С этим связано принятое в реализации модели допущение, что адсорбируются все пришедшие извне атомы.
При образовании зародыша, ему присваивается определенная ориентация оси С2, исходя из принципа минимизации энергии взаимодействия с окружением, как описано в пункте 1.
После адсорбции каждого атома происходит процесс диффузии. Учет возможности диффузии всех атомов, находящихся на поверхности пленки, сильно повышает вычислительную сложность задачи. Поэтому для диффузии выбираются не абсолютно случайные атомы с поверхности пленки висмута, а только те, что находятся в ближайшем окружении одного из нескольких недавно адсорбировавшихся атомов. При таком механизме «выборочной» диффузии учитываются последние изменения на поверхности пленки. Поэтому в модель введен коэффициент, определяющий, какая часть атомов поверхности будет диффундировать после адсорбции каждого атома. Изменяя значения этого коэффициента, можно добиться оптимального баланса между производительностью и точностью результатов модели.
3. Сравнение результатов моделирования с экспериментом
Исследование реальной структуры пленок висмута на подложках из слюды проводилось методом атомно-силовой микроскопии, выделение границ кристаллитов проводилось с помощью естественного оксидирования, взаимная ориентация кристаллитов определялась по фигурам роста [3]. Сравнение границ кристаллитов, полученных в модели, и границ кристаллитов реальных пленок показало высокую степень подобия (рис. 4). Сравнение закономерностей роста пленок висмута, полученных в модели и установленных экспериментально [5] позволяет определить их общность:
• С увеличением толщины растущей пленки наблюдается увеличение среднего размера кристаллитов, что происходит в процессе поглощения мелких кристаллитов более крупными.
• При увеличении температуры наблюдается увеличение размера кристаллитов, что связано в первую очередь с ускорением процесса диффузии.
• При увеличении скорости потока напыляемого материала наблюдается некоторое уменьшение среднего размера кристаллитов, связанное, скорее всего, с уменьшением относительной роли процессов диффузии.
• Текстура пленки, полученная с помощью компьютерной модели, полностью соответствует экспериментально наблюдаемой текстуре пленок висмута на слюде, полученных в широком интервале температур подложки <50-230° С).
• Важным результатом моделирования является рост пленок висмута на слюде с произвольными ориентациями кристаллографической оси С2 при превышении некоторого порогового значения температуры. Известно, что в процессе роста пленок методом термического напыления в вакууме такие кристаллиты появляются при температуре примерно 425 К [5].
Рис. 4а. Структура пленки, полученная с помощью компьютерной модели. Черным обозначены кристаллиты одной устойчивой ориентации, белым — другой
Рис. 4б. Реальная структура пленки Ы толщиной 300 нм. Кристаллиты одной ориентации подкрашены темным цветом, другой — светлым
4. Заключение
Предложенная модель процесса роста пленок висмута на подложках из слюды методом термического напыления позволяет оптимизировать процесс моделирования с учетом специфики физического процесса роста тонких пленок висмута. Исключение этапа реиспарения и оптимизация модели процесса диффузии позволили описать закономерности формирования и роста пленок висмута при уменьшении вычислительной сложности реализации модели.
Проведенное моделирование процесса роста пленок висмута на подложках из слюды, сравнение со структурой реальных пленок показало, что подложка вследствие своей структуры ориентирует первые слои зародышей висмута в наиболее энергетически выгодные антипараллельные ориентации, и эта структура пленки сохраняется по всей ее толщине.
Как следует из сравнения с экспериментом, модель, основанная на использовании принципа Больцмана и синусоидального закона изменения энергии границ раздела слюда — висмут и висмут — висмут в зависимости от кристаллографической ориентации, реализованная методом Монте-Карло, достаточно адекватно описывает явления, происходящие в процессе роста пленок висмута на слюде.
СПИСОК ЛИТЕРАТУРЫ
1. Александров Л. Н., Бочкова Р. В., Коган А. Н., Тихонова Н. Н. Моделирование роста и легирования полупроводниковых пленок методом Монте-Карло. Новосибирск: Наука, 1991. 168 с.
2. Бойко А. М., Сурис Р. А. Роль одномерной диффузии в модели роста поверхности кристалла Косселя // Физика и техника полупроводников. 2006. Т. 40. Вып. 3. С. 372-379.
3. Грабов В. М., Демидов Е. В., Комаров В. А., Климантов М. М. Атомно-силовая микроскопия декорированных оксидированием дефектов пленок висмута // Физика твердого тела. 2009. Т. 51. Вып. 4. С. 800-802.
4. Грабов В. М., Демидов Е. В., Комаров В. А. Атомно-силовая микроскопия пленок висмута // Физика твердого тела. 2008. Т. 50. Вып. 7. C. 1312-1316.
5. Грабов В. М., Демидов Е. В., Комаров В. А., Климантов М. М., Матвеев Д. Ю., Слепнев С. В.,
Усынин Е. В., Христич Е. Е., Константинов Е. В. Особенности структуры пленок висмута, полученных методом термического испарения в вакууме // Известия РГПУ им. А. И. Герцена. 2009. № 95. С. 105120.
6. Дубровский В. Г. Теория формирования эпитаксиальных наноструктур. М.: ФИЗМАТЛИТ, 2009. 352 с.
7. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1971. 328 с.
8. Киттель Ч. Введение в физику твердого тела. М.: Наука, 1978. 789 с.
9. Кукушкин С. А., Осипов А. В. Процессы конденсации тонких пленок // Успехи физических наук. 1998. Т. 168. № 10. С. 1083-1116.
10. Палатник Л. С., Фукс М. Я., Косевич В. М. Механизм образования и субструктура конденсированных пленок. М.: Наука, 1972. 320 с.
11. Рабинович В. А., Хавин З. Я. Краткий химический справочник Л.: Химия. 1977. 392 с.
REFERENCES
1. Aleksandrov L. N., Bochkova R. V., Kogan A. N., Tihonova N. N. Modelirovanie rosta i legirovanija poluprovodnikovyh plenok metodom Monte-Karlo. Novosibirsk: Nauka, 1991. 168 s.
2. Bojko A. M., Suris R. A. Rol' odnomernoj diffuzii v modeli rosta poverhnosti kristalla Kosselja // Fizika i tehnika poluprovodnikov. 2006. T. 40. Vyp. 3. S. 372-379.
3. Grabov V. M., Demidov E. V., Komarov V. A., KlimantovM. M. Atomno-silovaja mikroskopija dekori-rovannyh oksidirovaniem defektov plenok vismuta // Fizika tverdogo tela. 2009. T. 51. Vyp. 4. S. 800-802.
4. Grabov V. M., Demidov E. V., Komarov V. A. Atomno-silovaja mikroskopija plenok vismuta // Fizika tverdogo tela. 2008. T. 50. Vyp. 7. C. 1312-1316.
5. Grabov V. M., Demidov E. V., Komarov V. A., Klimantov M. M., Matveev D. JU., Slepnev S. V., Usynin E. V., Hristich E. E., Konstantinov E. V. Osobennosti struktury plenok vismuta, poluchennyh metodom ter-micheskogo isparenija v vakuume // Izvestija RGPU im. A. I. Gercena: Nauchnyj zhurnal. 2009. № 95. S. 105120.
6. Dubrovskij V. G Teorija formirovanija jepitaksial'nyh nanostruktur. M.: FIZMATLIT. 2009. 352 s.
7. Ermakov S. M. Metod Monte-Karlo i smezhnye voprosy. M.: Nauka, 1971. 328 s.
8. Kittel' CH. Vvedenie v fiziku tverdogo tela. M.: Nauka, 1978. 789 s.
9. Kukushkin S. A., Osipov A. V. Processy kondensacii tonkih plenok // Uspehi fizicheskih nauk. 1998. T. 168. № 10. S. 1083-1116.
10. Palatnik L. S., FuksM. JA., Kosevich V. M. Mehanizm obrazovanija i substruktura kondensirovannyh plenok. M.: Nauka, 1972. 320 s.
11. Rabinovich V. A., Havin Z. Ja. Kratkij himicheskij spravochnik L.: Himija, 1977. 392 s.
