Содержание к диссертации
Введение
Глава 1. Полосы роста в полупроводниковых кристаллах 14
1.1 Проблема полос роста в полупроводниковых кристаллах 15
1.2 Сегрегация примеси 16
1.3 Модели пограничного слоя 20
1.4 Модели, основанные на численном решении уравнений
гидродинамики 25
Выводы по главе 1 34
Глава 2. Математическая модель роста полупроводниковых кристаллов 36
2.1 Постановка задачи 38
2.2 Цилиндрические координаты 41
2.3 Характерные величины и числа подобия 44
2.4 Криволинейные координаты 46
2.5 Численная схема для уравнений Навье-Стокса в ортогональных цилиндрических координатах 51
2.5.1 Численная схема расщепления 53
2.5.2 Пространственные аппроксимации этапов расщепления
для уравнений движения 57
2.5.3 Пространственные аппроксимации этапов расщепления уравнения теплопроводности 68
2.5.4 Об аппроксимации и устойчивости численной схемы
2.6 Уравнение диффузии и задача Стефана 77
2.7 Возможное обобщение уравнений на трехмерный случай 78
Выводы по главе 2 81
Глава 3. Решение уравнений в криволинейных координатах 83
3.1 Преобразование уравнений к криволинейным координатам 84
3.1.1 Преобразование уравнений движения жидкости 86
3.1.1.1 Преобразование уравнения неразрывности 87
3.1.1.2 Преобразование производной по времени 88
3.1.1.3 Преобразование оператора переноса 90
3.1.1.4 Преобразование слагаемых градиента давления и
массовой силы 94
3.1.1.5 Преобразование слагаемых вязких сил 94
3.1.2 Окончательный вид уравнений в криволинейных координатах 98
3.2 Схема расщепления уравнений 101
3.2.1 Расщепление уравнений движения 102
3.2.2 Центральный этап и уравнение Пуассона 107
3.2.3 Расщепление уравнения теплопроводности 117
3.2.4 Этап «Z-теплопроводность» 119
3.2.5 Этап «R-теплопроводность» 122
3.2.6 Расщепление уравнения диффузии 124
3.3 Тестовая задача. Сравнение численного решения задачи
Стефана с аналитическим решением 127
3.3.1 Аналитическое решение задачи Стефана 127
3.3.2 Сравнение численного и аналитического решений 130
Выводы по главе 3 139
Глава 4. Моделирование нестационарных воздействий ростового оборудования 141
4.1 Аппаратура для выращивания кристаллов
методами вертикальной направленной кристаллиз ции 142
4.1.1 Контролирование температуры в печи 142
4.1.2 Неравномерное перемещение тигля 149
4.2 Моделирование нестационарных воздействий 151
4.2.1 Неравномерное перемещение тигля шаговым двигателем 154
4.2.2 Моделирование погрешности подержания температуры на нагревательных элементах 160
4.3 Сопоставление результатов с теоретическими моделями и экспериментальными данными 165
4.4 Требования, рекомендованные автором к ростовому оборудованию 167
Выводы по главе 4 170
Основные результаты работы 172
Список литературы. 174
- Модели пограничного слоя
- Численная схема для уравнений Навье-Стокса в ортогональных цилиндрических координатах
- Преобразование производной по времени
- Моделирование погрешности подержания температуры на нагревательных элементах
Модели пограничного слоя
Изделия и приборы, основанные на полупроводниковых технологиях, стали неотъемлемой частью современной жизни человека. Технические характеристики этих приборов во многом зависят от используемых для их производства материалов. Поэтому ключевой технологией среди всей последовательной цепочки производственного процесса создания микро- и оптоэлектронных устройств является выращивание монокристаллов полупроводников и полупроводниковых соединений.
Значительное место среди выращиваемых полупроводниковых кристаллов занимают монокристаллы кремния, используемые в производстве различных полупроводниковых приборов [13-15]. Кремниевые монокристаллы на сегодняшний день является основным материалом подложки при создании интегральных схем. Также широкое применение имеют некоторые соединения трех- и пятивалентных материалов A3B5, например арсенид галлия GaAs, обладающих полупроводниковыми свойствами. Областями применения приборов на основе таких соединений в основном являются оптические системы передачи данных и оптические устройства отображения информации.
Развитие микро- и оптоэлектронной техники невозможно без увеличения качества используемых полупроводниковых кристаллов и соединений. Поэтому получение качественных кристаллов необходимых размеров, чистоты и однородности электрофизических и оптических свойств является важной задачей. Для решения этой задачи разрабатываются новые и совершенствуются старые методы выращивания кристаллов. При этом мощным инструментом выступает математическое моделирование.
Свойства кристаллов в большей степени определяются распределением легирующих и остаточных примесей и дефектами структуры кристалла [16-17]. Для получаемых кристаллов характерны макроскопические (порядка 1 - 102 мм) и микроскопические (порядка 1 - 102 мкм) неоднородности примесей. В настоящее время существуют несколько технологических методик выращивания кристаллов, позволяющих влиять на макроскопические неоднородности посредством дополнительного перемешивания расплава: метод ACRT(Accelerated Crucible Rotation Technique) [18], AVT (Angular Vibration Technique) [19], выращивание кристаллов в магнитных полях [20], использование погружённого осциллятора [21,22] и др.
Существенной проблемой при выращивании полупроводниковых кристаллов и соединений является образование микроскопических концентрационных полос легирующей примеси (микросегрегаций). Из-за полос, наблюдаемых на микроскопических снимках поверхности кристалла после электрохимического травления (рис. 1.1), эти неоднородности принято называть полосами роста. Под полосами роста в настоящее время также подразумеваются и вариации состава сложных кристаллов или периодические флуктуации в структуре эвтектики. В результате, полосы роста приводят к неоднородности свойств получаемых подложек в пределах от одного микрометра до миллиметра [23, 24], что в дальнейшем оказывает негативное влияние на качество эпитаксиальных слоев и характеристики полупроводниковых приборов. ZnGeP2 , выращенного вертикальным методом Бриджмена В микро- и оптоэлектронике широко используются многокомпонентные полупроводниковые соединения или твердые растворы, такие как GaAsxPi-x, GaxInbxSb, CdxHgi_xTe и др. Концентрация компонентов варьируется и, тем самым, свойства получаемого твердого раствора оптимизируются в зависимости от условий применения. Так, меняя состав х твердого раствора CdxHgi_xTe в пределах 0,2 - 0,5, можно получать материал, пригодный для детектирования ИК-диапазона в трёх диапазонах длины волны: 1,5 - 2,5; 3 - 5 и 8 - 14 мкм [25]. При этом современные мировые требования к однородности состава являются достаточно высокими, см. таблицу 1.1 [9]. Например, для полупроводниковых подложек из CdxHgbxTe необходимо, чтобы максимально допустимая неоднородность состава составляла AJC/JC«10 4-10 5. Образующиеся в кристалле полосы роста часто нарушают это требование, что приводит к уменьшению эффективности полупроводниковых приборов или полной непригодности выращенных кристаллов.
Численная схема для уравнений Навье-Стокса в ортогональных цилиндрических координатах
В системе уравнений (2.1.1) - (2.1.5) приняты следующие обозначения: v - вектор скорости, р - давление, Т - температура, р -плотность, р0 - средняя плотность, g - ускорение свободного падения, JU -коэффициент динамической вязкости, Ль - коэффициент температуропроводности расплава, D - коэффициент диффузии, С - концентрация примеси. Уравнения записаны в предположении, что XL и D являются постоянными величинами.
Сущность приближения Буссинеска заключается в учете сжимаемости среды только в массовой силе уравнения (2.1.1), а во всех других уравнениях плотность среды полагается постоянной. Несмотря на кажущуюся непоследовательность, это приближение очень хорошо себя зарекомендовало, и большая часть вычислительной информации об естественной конвекции в жидкостях получена с помощью решения уравнений в приближении Буссинеска.
В кристалле и тигле процесс переноса тепла описывается дифференциальным уравнением теплопроводности СР g где т, т - температуры на фронте кристаллизации в кристалле и в расплаве, Гп температура плавления, у - удельная теплота кристаллизации, СР - удельная теплоемкость, vg - скорость роста кристалла. Первое условие Стефана (2.1.8) означает непрерывность температуры на фронте кристаллизации и её равенство температуре плавления. Второе условие (2.1.9) означает, что разность тепловых потоков на фронте равна теплоте кристаллизации. Это условие записано в предположении, что плотности кристалла и расплава одинаковы.
При расчете распределения примеси на фронте кристаллизации учитывается эффект сегрегации. Дело в том, что концентрация примеси в закристаллизовавшейся части отличается от концентрации примеси, которая была в момент кристаллизации в расплаве. Отношение этих концентраций в равновесном состоянии кристалла и расплава и при малых скоростях роста называют коэффициентом сегрегации примеси = —, где Cv- концентрация в
Выращивание полупроводниковых кристаллов из расплава методами вертикальной направленной кристаллизации происходит в цилиндрических тиглях. Таким образом, область решения модельной задачи также представляет собой цилиндр, и все приведённый выше уравнения следует записать в цилиндрической системе координат.
Для перехода к цилиндрической системе координат необходимо определить градиент и Наряду с граничными условиями на внутренней подвижной границе раздела фаз (2.1.8 - 2.1.10) в задаче используются следующие граничные условия. Для вектора скорости на границе с тиглем и фронтом кристаллизации ставится условие непротекания и прилипания: и = 0, w = 0; на оси г = 0 - равенство нулю функции «вихрь»: и = 0, — = 0. На нагревателях r = R+H задается температурный режим T = THm(z), на верхнем и нижнем торцах - Т = Тгор(г) и Т = Твоггом(г) соответственно. На границе кристалла и расплава с тиглем r = R ставится условие равенства радиальных тепловых потоков, а на оси г = 0 - условие симметрии — = 0. Для концентрации примеси на границе с тиглем и на верхней границе расплава ставится условие отсутствие диффузионного потока, а на оси = 0. 2.3 Характерные величины и числа подобия Введём безразмерные переменные, определяемые по следующему правилу: r = Zr , z = Zz , где Z - характерный размер области, v = Uv , t = —f. В качестве характерной скорости выбрано U = Jj3-AT-g-Z, где р - коэффициент температурного расширения материала, AT - перепад температуры в области. Для плотности, температуры, концентрации примеси и давления получим Р = Р0р , Т = АТ-Г + Т0, С = АС-С + С0, р = Р-р + Р0. Здесь Р0- характерная плотность, АС - характерный перепад концентрации примеси в области, Г0 и С0 -размерные значения температуры и концентрации, которые равны температуре кристаллизации и концентрации примеси соответственно, P = p0U , гидростатическое давление Р0 =p0gz.
При решении двухфазной задачи Стефана фронт кристаллизации является криволинейным и подвижным, что приводит к необходимости использования криволинейных координат. Введём функции преобразования координат от переменных (r,z,t) к переменным (,/7, г) следующим образом Преобразование (2.4.1) таково, что область решения в плоскости {r,z) перейдет в единичный квадрат в плоскости (,77), а линия фронта кристаллизации перейдёт в линию r, = rf=comt. Таким образом, вместо неравномерной, криволинейной, нестационарной сетки в физической области мы получаем равномерную, прямоугольную, стационарную сетку в вычислительной области (рис. 2.3). Сетка по г выбрана равномерная и зависит только от координаты , а для координаты z сетка зависит от , г] и т. Узлы сетки по z сгущаются в расплаве вблизи фронта кристаллизации z = q (r,i) и верхней границы z = \. В твердом теле кристалла сетка по переменной z равномерная для каждого г (или, соответственно, ).
Преобразование производной по времени
В работе [3] Фавье и др. представили сравнение своей модели с моделями Ван Рана [2] и Вильсон [29-32]. Сравнение моделей проводилось при параметрах, соответствующих сегрегации легирующих элементов кислорода и углерода в кремнии. Результаты расчета по разным моделям приведены на рисунке 1.4. На рисунке представлены относительные концентрации CS/CL в зависимости от координаты в направлении роста для части кристалла, выращенной за один цикл изменения скорости. Результаты для стационарной модели БПС оцениваются от точки к точке. По рисунку видно, что стационарная модель дает слишком большие амплитуды флуктуаций и поэтому не годится для количественного сопоставления. Однако три другие нестационарные модели по величине микронеоднородности прекрасно согласуются между собой. Рисунок 1.4. Относительные концентрации CS/CL кислорода (а) и углерода (б) в кремнии в зависимости от координаты в направлении роста для части кристалла, выращенной за один цикл изменения скорости. Расчёт по различным моделям
Г. Мюллером в монографии [26] было выполнено сопоставление результатов теоретической модели Ван Рана с экспериментальными данными, доступными из работы [28]. Для сравнения из работы [28] были взяты два экспериментальных результата: микроскопические изменения концентрации кислорода в легированном бором кристалле кремния, выращенном методом Чохральского (кристалл совершал 4 оборота в минуту, скорость вытягивания R0 составляла 40 мм/час, рис. 1.5) и микроскопические изменения скорости роста для эксперимента, где скорость вытягивания R0 составляла 40 мм/час, а кристалл совершал 5 оборотов в минуту (рис. 1.6). По колебаниям скорости роста на рисунке 1.6 был оценен параметр А = 4. Величина К0Г1/2 была известна из параметров эксперимента и составила 0.4 мм/мин. По результатам, полученным Ван Раном и представленным на рисунке 1.3, была вычислена величина микроскопической неоднородности, которая составила 10-15%, что согласуется с экспериментальными результатами, приведенными на рисунке 1.5. Таким образом, было сделано заключение, что указанные модели можно использовать для описания микросегрегаций, вызванных флуктуациями скорости роста. Действительно, если известна информация о колебаниях скорости роста, данные модели возможно применить. Однако откуда эту информацию взять? Для реальных экспериментов такая информация, как правило, не доступна. Поэтому использование моделей пограничного слоя ограничивается необходимостью расчета скорости роста с помощь других математических моделей. Таким образом, возможность расчёта скорости роста кристалла оказывается ключевым фактором для моделирования микросегрегаций.
Ранее в разделе 1.2 Главы 1 было отмечено, что колебания скорости роста могут происходить по различным причинам, в первую очередь в результате нестационарных воздействий со стороны оборудования или в результате нестационарной конвекции в расплаве. Поэтому математические модели, которые способны учитывать такие нестационарные воздействия, представляют большой интерес для расчёта реальной величины микроскопической неоднородности.
На сегодняшний день для получения полупроводниковых кристаллов широко используются два метода выращивания кристаллов: метод Чохральского (и его модификации) и методы вертикальной направленной кристаллизации.
Соответственно, в зависимости от материала и требований может применяться тот или иной метод роста. Например, для выращивания кристаллов кремния в настоящее время используется только метод Чохральского. Он экономически выгоден и позволяет получать кристаллы большого диаметра и с низкой плотностью дислокаций. Но для получения кристаллов полупроводниковых соединений (таких, как арсенид галлия или фосфид индия) он не всегда оказывается удачным, так как получаемые кристаллы имею высокую плотность дислокаций. Для примера рассмотрим монокристаллы арсенида галлия GaAs. Более 95% рынка GaAs составляют кристаллы двух типов: - полуизолирующий GaAs с удельным сопротивлением более 107 Ом, используемый при производстве высокочастотных ИС (до 24 ГГц) и дискретных приборов. Такие кристаллы, как правило, получают методом Чохральского; - Сильно легированный кремнием Si (1017–1018 см-3) GaAs с низкой плотностью дислокаций, применяемый при изготовлении светодиодов и лазеров. Такие монокристаллы GaAs должны обладать низкой плотностью дислокаций и, поэтому, для его производства используются методы вертикальной направленной кристаллизации.
При выращивании кристаллов методом Чохральского главной причиной колебаний скорости роста является интенсивная нестационарная конвекция в расплаве [26]. Одновременный расчёт трехмерной нестационарной конвекции и скорости роста кристалла представляет собой очень сложную задачу. На сегодняшний день не существует математических моделей, способных рассчитать микросегрегации для метода Чохральского. При выращивании кристаллов методом вертикальной направленной кристаллизации с затравкой внизу, напротив, конвекция обычно имеет стационарный ламинарный характер. Поэтому полос, связанных с нестационарной конвекцией, здесь нет. Однако, как было показано в экспериментальной работе [39], колебания скорости роста кристалла и микросегрегации все же возникают, см. рисунок 1.7.
Авторы работы утверждают, что основной причиной осцилляций скорости роста являются флуктуации температуры, вызванные погрешностями в системе регулирования технологической установки. В дополнение заметим, что в научной литературе можно встретить и другие работы, где также подтверждается этот факт, например [8, 26, 37]. Поэтому в методах вертикальной направленной кристаллизации основной причиной колебаний скорости рост кристалла оказываются именно нестационарные воздействия в оборудовании. Здесь нужно особенно выделить такие воздействия как неравномерное перемещение тигля относительно нагревательных элементов шаговым двигателем и погрешность в поддержании требуемой температуры в печи.
Моделирование погрешности подержания температуры на нагревательных элементах
В процессе направленной кристаллизации вещество внутри тигля находится в двух (жидком и твердом) состояниях, между которыми существует четкая граница раздела фаз. В расплаве тепло- и массоперенос осуществляется естественной конвекцией, теплопроводностью и диффузией. Конвекция, в свою очередь, образуется в результате градиентов температуры и концентрации в поле силы тяжести. На фронте кристаллизации происходит выделение скрытой теплоты кристаллизации. В кристалле тепло- и массоперенос осуществляется теплопроводностью и диффузией соответственно. Однако следует отметить, что коэффициент диффузии в кристалле примерно на 4 порядка меньше чем в расплаве, поэтому диффузией в кристалле обычно пренебрегают. Таким образом, задача моделирования процесса направленной кристаллизации заключается в совместном решении нестационарной задачи Стефана и системы уравнений гидродинамики в расплаве или системы уравнений Навье-Стокса. При этом систему уравнений Навье-Стокса, как правило, рассматривают в приближении Буссинеска [8,37, 52-54].
Для решения такой задачи, в диссертационной работе используется численная схема, ранее разработанная научным руководителем В.А. Гончаровым. Многочисленные аспекты этой схемы были опубликованы в работах [10] и [11]. Однако для целостной формулировки разностной схемы, которая используется в данной диссертационной работе, необходимо сделать некоторые повторения, для того чтобы схема была описана полностью.
Для численного решения уравнений Навье-Стокса в приближении Буссинеска существуют два основных подхода, связанных с выбором расчётных переменных. При решении уравнении в естественных переменных скорость давление основным преимуществом является простота постановки граничных условий для поля скорости. Существенной трудностью при этом является необходимость решать уравнение Пуассона с условиями Неймана. Для другого подхода, использующего переменные вихрь-функция тока, указанной вычислительной трудности удается избежать. Для функции тока также необходимо решать уравнение Пуассона, но уже с граничными условиями Дирихле, что является значительно более легкой вычислительной задачей. В то же время такой подход имеет существенный недостаток, связанный с появлением дополнительного граничного условия для функции вихря, что приводит к дополнительным погрешностям в решении. Разработанный В.А. Гончаровым и подробно описанный в работе [10] численный метод, основанный на специальном расщеплении и использовании разнесённых пространственных сток, позволил совместить удачные стороны обоих способов расчёта конвекции. В результате, используются естественные переменные, однако вместо уравнения Пуассона с условиями Неймана решается уравнение Пуассона с условиями Дирихле для разностного аналога функции тока.
Консервативность численной схемы и использование неравномерных пространственных сеток позволяют получать достаточно точное решение задач конвекции на сравнительно небольшом числе узлов. Это подтверждается решением задачи [55], которая является общепринятым тестом для двумерных численных схем в этой области.
Специфика численного метода, которая заключается в специальном расщеплении по физическим процессам, а там, где это оправдано, и по пространственным переменным, позволила разработать принципиально новый способ решения задачи Стефана. Благодаря сквозному счёту уравнения теплопроводности в кристалле и расплаве удается на каждом шаге по времени вычислять изменение температурного поля в этих областях и мгновенную скорость роста кристалла [11]. Обратим внимание, что для рассматриваемой проблемы расчёта величины микроскопической неоднородности это является ключевым фактором. Следует отметить, что способ расчёта температурного поля и скорости роста является квазиодномерным и может использоваться для решения как двухмерных, так и трехмерных задач.
Для моделирования процесса роста кристалла выделим области расплава, кристалла и тигля (рис. 2.2). В области расплава учитывается течение жидкости под влиянием естественной конвекции. В кристалле и тигле рассчитывается перенос тепла в результате теплопроводности материалов. Перенос легирующей примеси рассматривается только в расплаве в результате диффузии и конвективных течений. На фронте кристаллизации ставится задача о фазовом переходе, или задача Стефана, и учитывается эффект сегрегации примеси. Рассматривается рост кристаллов с малым содержанием примеси, поэтому полагается, что температура фазовой границы не зависит от концентрации примеси. Также не учитывается зависимость плотности от концентрации примеси. Задача рассматривается в двумерной осесимметричной постановке.
