УДК 548.4 + 519.24
С. С. Гуськов, М. А. Фаддеев, Е. В. Чупрунов
МОДЕЛИРОВАНИЕ РОСТА ДВУХКОМПОНЕНТНОГО КРИСТАЛЛА СО СЛОИСТЫМ РАСПРЕДЕЛЕНИЕМ ПРИМЕСИ
Предлагается оригинальный алгоритм стохастического моделирования роста двухкомпонентного кристалла, в состав которого может входить как основное вещество, так и примесь. Показано, что в зависимости от ряда параметров модели возможно получение как равномерного, так и квазипериодического распределения примеси, т.е. формирование зонарной микроструктуры кристалла. Проведен анализ условий, приводящих к неравномерному распределению примеси. Получено качественное согласие результатов моделирования с экспериментальными данными.
Введение
Большинство природных и синтетических монокристаллов содержат дефекты, известные как зонарная микроструктура, полосчатость или слоистость. Эти дефекты проявляются в виде чередующихся слоев неодинакового качества, параллельных фронту кристаллизации. Характерный период зонарной структуры составляет единицы и десятки микрометров. Согласно современным представлениям [1], различие физических характеристик отдельных слоев обусловлено неравномерным распределением примеси. Для формирования полосчатости достаточно малой концентрации примеси (—0,01%). Зо-нарная структура существенно влияет на оптические свойства кристаллов, что необходимо учитывать, в частности, при использовании их в лазерной технике в качестве разного рода оптических элементов.
Целью данной работы является исследование зонарной структуры кристаллов (ЗСК) и анализ факторов, влияющих на зонарность. Предлагается механизм образования ЗСК за счет накопления примеси в среде у поверхности растущей грани и последующего лавинообразного встраивания примеси при превышении некоторого критического уровня ее концентрации. Разработана математическая модель роста двухкомпонентного кристалла. С помощью численного статистического моделирования показана возможность формирования как практически однородного распределения примеси в кристалле, так и квазипериодической зонарной микроструктуры.
1. Зонарная микроструктура кристаллов
Зонарная микроструктура - широко распространенное явление [1, 2], которое характерно для многих естественных минералов (рис. 1,а), а также для кристаллов, выращенных из расплава (рис. 1,б) и раствора (рис. 1,в). Тем не менее причины и механизмы неравномерного распределения примеси изучены слабо, общий метод описания зонарной структуры не разработан. Попытки адекватного объяснения ЗСК обычно основаны на одном из двух различных предположений: 1) зонарное распределение - результат внутренних процессов роста кристалла [3]; 2) зонарное распределение отражает изменения во внешней среде в течение роста кристалла [1, 4]. Первая гипотеза подразумевает, что ЗСК генерируется без влияния временных зависимостей внешних условий, т.е. путем самоорганизации. Второе предположение исходит из связи ЗСК с колебанием внешних условий роста (температуры, пересыщения) за счет конвекционных потоков.
а) б) в)
Рис. 1 Зонарная микроструктура кристаллов: а - гранат, увеличение х10 [2]; б - Ти8Ъ, увеличение х95 [4]; в - КОР, увеличение х160 [5]
ЗСК содержит ценную информацию о процессах роста кристалла и изменениях в среде кристаллизации. Исследование ЗСК позволяет оценивать качество ряда кристаллов, важных в лазерном приборостроении [5], и в нужном направлении совершенствовать технологию их выращивания. Наличие микрослоистой структуры может приводить и к возникновению новых свойств кристалла [6].
2. Математическая модель роста кристалла с примесью
Процессы, происходящие при росте кристаллов, очень сложны и могут быть описаны функциями, зависящими от многих параметров, характеризующих внешние условия. По этой причине одним из наиболее перспективных способов изучения этих процессов является метод численного статистического моделирования (метод Монте-Карло) [7].
Будем рассматривать кристаллизацию как переход ростовых единиц (РЕ) - атомов, молекул и т.п. - из маточной среды в кристалл. В данной работе моделируется рост грани (001) кристалла с простой кубической решеткой (кристалл Косселя) с использованием циклических граничных условий; при этом рассматривается 2 типа РЕ (основное вещество и примесь). РЕ обоих типов могут встраиваться только в узлы решетки; используется приближение ближайших соседей - рассматривается 6 соседей каждой РЕ.
Пусть узел на поверхности кристалла имеет координаты (х, у, г). Возможны 2 типа элементарных событий (рис. 2): присоединение РЕ определенного типа (г = 1, 2) в узел (х, у, г+1) и отрыв от поверхности кристалла РЕ, находящейся в узле (х, у, г). Множество W возможных элементарных событий включает 3 элемента:
W = ^1+, W2+, W-}. (1)
Пусть потенциальная энергия РЕ в кристалле #1, а в маточной среде -Е2 (очевидно, что при росте кристалла #1 < Е2). Вероятности нахождения РЕ в соответствующих средах определяются соотношениями:
'2
р1 = А ехРІ- кТ
, Р2 = Аехр| -2 1 кТ
(2)
где к - постоянная Больцмана; Т - абсолютная температура; А - нормировочный коэффициент.
Рис. 2 Присоединение (1) и отрыв (2) РЕ от поверхности кристалла
Пусть Pi и Pj - вероятности нахождения РЕ в позициях I и у. Тогда вероятность перехода РЕ из положения i в положение у равна
У = В-
(3)
где В - нормировочная постоянная.
Если г соответствует кристаллу, а у - среде, то переход у ^ г представляет собой присоединение РЕ к кристаллу. Обратный переход соответствует десорбции (отрыву) РЕ. Подставим (2) в (3) и обозначим
АЕ = Е2 — Е > 0 .
(4)
Тогда вероятности присоединения и отрыва выразятся следующими формулами:
W+= В----------------^-г
1 + ехр(—АЕ / кТ)
, W— = В
1
1 + ех р(АЕ / кТ)'
(5)
Способ расчета энергии РЕ в среде и кристалле определяется используемыми моделями среды и кристалла. В нашей модели энергия данной РЕ в среде полагается постоянной, а в кристалле - зависит от числа и типа соседей:
Е1 = —1е, Е2 = Е0,
(6)
здесь I - число соседей данной РЕ; е > 0 - энергия связи двух РЕ.
Изменение величины энергии #0 позволяет управлять скоростью кристаллизации. Следовательно, энергию #0 можно представить как некоторую функцию параметра о, который описывает отклонение системы «кристалл -среда» от равновесия. Параметр о определяется в случае роста из паровой фазы давлением пара, в случае роста из раствора - концентрацией кристал-
лизующегося вещества в растворе, при росте из расплава - температурой расплава [5].
Путем несложных преобразований выражений (5) с учетом (6) получим
W+= В----------------------Ц-г, W—= В--------1---- ------, (7)
1 + / (а)ехр(- 1е / кТ) 1 + (/ (Ст))—1 ехр(е / кТ)
где введено обозначение
/ (а) = ехр(—Е0(а)/ кТ). (8)
Для анализа общих свойств зависимости #0(о) воспользуемся принципом детального равновесия: при нулевом отклонении системы «кристалл -среда» от равновесия кристалл расти не должен. Это означает, что при о = 0 вероятности W+ и W- равны. Тогда из (4), (5) и (6) следует, что при о = 0 энергия #0 = - <1>е. Так как рассматривается кристалл Косселя и используется приближение ближайших соседей, то среднее число соседей частицы <1> = 3 и энергия #0 = - 3е. При возрастании параметра о значение #0 должно увеличиваться, что соответствует ускорению роста кристалла. При о < 0 кристалл растворяется. Представим зависимость #0(о) в следующем виде:
Е0(а) = 3е(Ш(а) — 1). (9)
В случае двухкомпонентного кристалла (2 типа РЕ) вероятность присоединения должна зависеть от концентрации данного типа РЕ в слое среды у поверхности кристалла. Влиянием концентрации примеси на вероятность отрыва пренебрежем:
ВС В
^ = 1+„, Е /.т), W—г - W—=-----------------------1--------------. (10)
1 + /(а)ехр(Е1г-/ кТ) 1 + (/(а)) 1 ехр(- Е1г /кТ)
Здесь индекс г обозначает тип РЕ, г = 1 для основного вещества, г = 2 для примеси. Для двухкомпонентного кристалла энергию Е1г- РЕ г'-го типа в кристалле представим, аналогично (6), следующим соотношением:
Е1г =—11ег1 — 12ег 2, (11)
где 1к - число соседей типа к; егк > 0 - энергия связи РЕ г'-го и к-го типа. Для кристалла с примесью в выражении (9) для #0(о) энергия е заменяется на энергию связи РЕ основного вещества ец.
Согласно (10), с увеличением о повышаются вероятности присоединения W+i и понижается вероятность W- отрыва РЕ, т.е. увеличивается скорость роста. С увеличением разности энергий РЕ в среде и в кристалле вероятность присоединения растет, вероятность отрыва снижается. С ростом температуры Т понижается вероятность присоединения и увеличивается вероятность отрыва РЕ, т.е. уменьшается скорость роста, что соответствует экспериментальным данным.
Базовый алгоритм моделирования заключается в повторении следующих операций:
1. Выбор текущего узла на поверхности кристалла.
2. Расчет нормированных вероятностей возможных событий (10).
3. Выбор элементарного события в зависимости от значения псевдослучайного числа.
4. Выполнение события.
При использовании численного стохастического моделирования время г процессов принято измерять в шагах алгоритма Монте-Карло (МКШ). Один МКШ реализуется при переборе всех узлов решетки на поверхности рассматриваемого участка грани. Расчет продолжается в течение заданного числа МКШ.
Рассмотрим слой маточной среды у поверхности кристалла - приповерхностный слой (ПС). Примем суммарную концентрацию РЕ рассматриваемых двух типов в ПС за единицу. Обозначим концентрацию примеси С, тогда концентрация основного вещества равна 1 - С. Перед началом расчетов
задается средняя концентрация примеси в среде С0. После заданного числа МКШ производится пересчет концентрации основного вещества и примеси в ПС. Пусть ПС содержит фиксированное число щ РЕ, причем из них щ типа г. Тогда концентрация примеси определяется соотношением
С = Щ2. (12)
П0
Пусть за время Дг к кристаллу присоединилось к0 РЕ, причем среди них кг РЕ типа г. Присоединившиеся РЕ покинули ПС, на смену им из среды кристаллизации в ПС поступило к0С0 РЕ примеси и £ф(1-С0) РЕ основного вещества. Таким образом, новая концентрация в ПС определяется через старую следующим образом:
С(7+1) = п2 — к2 + к0С0 = с(7) + к0С0 — к2 , с(0) = С(0, (13)
п0 п0 ’ ’
где 7 - порядковый номер этапа алгоритма.
Для проверки работы модели после выбора оптимального времени расчета и размеров области расчета было проведено исследование роста при разных входных параметрах. В частности рассчитывался коэффициент присоединения РЕ, который характеризует скорость роста грани и определяется как отношение общего числа присоединившихся РЕ п к полному числу итераций т: V = п/т. Зависимости коэффициентов присоединения от отклонения от равновесия о при разных значениях ец/кТ для однокомпонентного (без примеси) кристалла приведены на рисунке 3.
Для уменьшения дисперсии результатов производилось трехкратное повторение каждого компьютерного эксперимента с последующим усреднением характеристик. Установлено, что результаты моделирования при различных параметрах кристалла качественно совпадают с экспериментальными данными. На каждой кривой v(о) наблюдаются практически ли-
нейные участки, характерные для роста шероховатых граней реальных кристаллов [5].
V
1.0
0.8
0.6
0.4
0.2
0.0
0.0 0.2 0.4 0.6 0.8 1.0
Рис. 3 Зависимости у(с) при разных е11/кТ (С0 = 0)
3. Распределение примеси в растущем кристалле
Для анализа условий неравномерного распределения примеси проведено моделирование роста двухкомпонентного кристалла. Энергии связи в твердых телах принимают значения от сотых (молекулярные кристаллы) до единиц (ионные кристаллы) электрон-вольт [8]. Следовательно, реальным
кристаллам при Т = 300 К соответствует диапазон 0,5 < ец/кТ < 500. Отметим, что при разработке теории дислокационного роста для оценок использовалась величина е/кТ = 4 [9]. Энергия связи РЕ основного вещества с РЕ примеси е]_2 и энергия связи РЕ примеси е22 не должны отличаться от ец более чем на два порядка: 0,01 < е^/ец < 100, 0,01 < е22/ец < 100.
Зафиксируем значения о, С0 и ец/кТ: о = 0,2, С0 = 0,01, ец/кТ = 4, и рассмотрим распределения примеси при разных отношениях е12/е11 и е22/е11. Возможна реализация четырех вариантов поведения системы.
Варианту 1 (см. рис. 4,а) соответствует равномерное распределение примеси в кристалле без накопления и без дефицита примеси в среде у фронта кристаллизации. Варианты 2 и 3 - это равномерное распределение примеси в кристалле с накоплением и с дефицитом примеси соответственно в ПС среды.
Вариант 4 (рис. 4,б) представляет собой квазипериодическое распределение примеси в кристалле с соответствующими колебаниями концентрации примеси в среде. Механизм образования слоев примеси следующий: примесь в среде у поверхности растущей грани накапливается и при превышении критического уровня концентрации лавинообразно присоединяется к кристаллу.
Изменяя е 12/Е 11 и £22/£11 при нескольких постоянных значениях ец/кТ, построим ряд энергетических диаграмм. При этом четыре рассмотренных варианта распределения примеси будут наглядно представлены в виде областей
в осях е 12/е 11 - Е22/Е11 (рис. 5,а).
0.05 т С 0.04 0,03 0,02 0,01
0
0,251 С 0.20 -0,150,100.05-
0,05 і С 0,04 -0,03 -0,02 -0,01
г, а0
а)
0
500
1000
1500
о
500
1000
1500
б)
Рис. 4 Распределение примеси: а - вариант 1; б - вариант 4. Слева - распределение примеси в кристалле, справа - зависимость концентрации примеси в ПС от времени; а0 - параметр элементарной ячейки
г, МКШ
2000
В зависимости от величины єц/кТ область квазипериодического распределения примеси (область 4) меняет конфигурацию и размеры, а при єц/кТ < 1 вообще исчезает с энергетической диаграммы.
Размер и положение областей с разными вариантами распределения примеси довольно чувствительны к изменению отклонения от равновесия о
(рис. 5,б) и средней концентрации примеси в маточной среде Со (рис. 5,в).
В результате моделирования выявлены следующие закономерности поведения параметров слоистого распределения примеси. Во-первых, при увеличении параметра о уменьшаются период и толщина легированных примесью слоев. Во-вторых, снижение концентрации примеси в среде Со ведет к уменьшению содержания примеси в слоях и к увеличению их периода, т.к. примесь медленнее накапливается в ПС, и пока ее концентрация увеличивается до критической, успевает вырасти более толстый чистый слой. Таким образом, при изменении параметров маточной среды возможно возникновение квазипериодического распределения примеси с разными количественными характеристиками.
е„/кТ=2
Єі1/кТ=4
3
■1
2
4
0.01 01 об 1 г
е^/е 11
100
10
5
2
1
0.5
0.1
0.01 10 100
3
1
I I
2 і— / \ ■
4 / 1
100
10
5
2 *-
Ф
1 *- Ф
0.5
0.1
0.01
0.01 0.1 0.5 1 2
а) Єц/Єн
10 100
и = 0.1
сг= 0.3
3
1
2
і ■
0.01 0.1 0.5 1 2
е22^ 11
2
1
г-
0.5
и
0.1
0.01
10 100
б)
3
1
2
4
0.01 0.1 0.5 1 2
егг/е 11
2
1
0.5
ф
0.1
0.01
10 100
С0= 0.002
3
1
2
4
0.01 0.1 0.5 1 2
Є22/Є 11
2
1
*-
О.Ь
ф
0.1
0.01
10 100
С = 0.05
1
/
2 I
і ■ 3
2
1
*-
0.5
ф
0.1
0.01
в)
0.01 0.1 0.5 1 2
егг/е 11
10 100
Рис. 5 Области с разным характером распределения примеси: а - при различных значениях ец/кТ; б - при различных значениях о; в - при различных значениях С0. Цифровые обозначения областей соответствуют номерам вариантов распределения примеси
Заключение
Разработанная модель роста двухкомпонентного кристалла позволяет объяснить формирование квазипериодической зонарной микроструктуры, исходя из предположения о возможности накопления примеси перед фронтом кристаллизации. При этом образуется кристалл, в котором высокосовершенные слои основного вещества разделяются более тонкими прослойками с
повышенным содержанием примеси. Показана возможность применения модели для исследования условий возникновения зонарности и расчета характеристик зонарного распределения примеси.
Список литературы
1. Современная кристаллография / под ред. А. А. Чернов [и др.]. - М. : Наука, 1980. -Т. 3 : Образование кристаллов. - 408 с.
2. Holten, T. Statistical characteristics and origin of oscillatory zoning in crystals / T. Holten, B. Jamtveit, P. Meakin [et al.] // American Mineralogist. - 1997. - V. 82. -P. 596-606.
3. Smolsky, I. L. X-ray topographic study of striation formation in layer growth of crystals from solutions / I. L. Smolsky, A. E. Voloshin, N. P. Zaitseva [et al.] // Philosophical Transactions of the Royal Society of London. - 1999. - V. 357. - P. 2631-2649.
4. Проблемы роста кристаллов / под ред. Н. Н. Шефталь и Е. И. Гиваргизова. - М. : Мир, 1968. - 392 с.
5. Портнов, В. Н. Возникновение и рост кристаллов / В. Н. Портнов, Е. В. Чупрунов. - М. : Физматлит, 2006. - 328 с.
6. Зайцева, Е. В. Рентгеновская дифракция при воздействии лазерного излучения на кристаллы эпсомита, выращенные в присутствии буры / Е. В. Зайцева, В. Н. Портнов, М. А. Фаддеев, Е. В. Чупрунов // Кристаллография. - 1997. - № 6. -42 т. - С. 969-971.
7. Методы Монте-Карло в статистической физике / под ред. К. Биндер - М. : Мир, 1982. - 444 с.
8. Павлов, П. В. Физика твердого тела / П. В. Павлов, А. Ф. Хохлов. - Н. Новгород : ННГУ, 1993. - 490 с.
9. Элементарные процессы роста кристаллов / ред. Г. Г. Леммлейн, А. А. Чернов. -М. : Издательство иностранной литературы, 1959. - 300 с.