Известия РАН. Механика твердого тела, 2020, № 5, стр. 51-63

1. Введение. Полупроводниковый кремний и основная технология его выращивания методом направленной кристаллизации из расплава по Чохральскому являются основой широкого развития твердотельной электроники. Компоненты установки для выращивания монокристаллов кремния методом Чохральского имеют осевую симметрию. Схема установки для выращивания показана на рис. 1: 1 – кристалл, 2 – расплав, 3 – фронт кристаллизации (ФК), 4 – тигель, 5 – подставка, 6 – нагреватель, 7 – теплоизоляция, 8 – тепловой экран, 9 – водоохлаждаемый корпус; габаритные размеры: H = 146 см и D = 86 см. Цилиндрические слитки монокристаллического кремния длиной L_к = 100–150 см и диаметром 15 и 20 см вытягиваются со скоростью V снизу вверх из расплава, содержащегося в графитовом тигле, окруженном резистивным высокотемпературным нагревателем, при этом вся конструкция помещена в водоохлаждаемую камеру и дополнена тепловыми экранами для управления температурой в растущем монокристалле. Кристалл вращается с угловой скоростью Ω_k, а тигель с подставкой вращается в противоположную сторону со скоростью Ω_Т.

Рис. 1

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

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

Применение современных моделей и методов моделирования положено в основу рассматриваемого направления, развиваемого в ИПМех РАН. Создаваемые механико-математические модели верифицируются по экспериментальным данным и используются для оптимизации конструкций тепловых узлов большегрузных промышленных ростовых установок и температурно-динамических условий выращивания монокристаллов и температурно-временных режимов термообработки пластин, гарантирующих контроль природы, содержания, размеров и характера распределения в них микродефектов [1–3].

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

Учет конвекции в расплаве осуществляется на основе сопряжения радиационно-кондуктивной модели и модели кондуктивно-конвективных процессов в системе кристалл–расплав с учетом выделения скрытой теплоты кристаллизации Q_H. Для этого на границах кристалла и расплава задаются распределения температуры, полученные в результате расчета по радиационно-кондуктивной модели для всего объема теплового узла.

Гидродинамическая модель строится на основе численного решения полных трехмерных уравнений Навье–Стокса совместно с уравнениями конвективного переноса тепла в расплаве за счет вращения кристалла и тигля, а также тепловой гравитационной конвекции. Параметры гидродинамической модели записываются в безразмерном виде: Re = Ω_kR_k/ν – число Рейнольдса, Ros = Ω_k/Ω_Т – число Россби, Gr = gβ_T$R_{k}^{3}$ΔT/ν² – число Грасгофа, Pr = ν/a – число Прандтля. Здесь: R_k – радиус кристалла, Ω_k, Ω_Т – угловые скорости вращения кристалла и тигля, ν – кинематическая вязкость расплава, g – ускорение силы тяжести, β_Т = –ρ^–1∂ρ/∂Т – коэффициент теплового расширения, ρ – плотность, ΔТ – максимум отклонения температуры T от температуры кристаллизации Т_k.

Важную роль в адекватности такого сопряженного расчета играют теплофизические параметры твердотельных компонент (табл. 1) и расплава кремния (табл. 2), приведенные с указанием номера соответствующей компоненты теплового узла на рис. 1.

Рассмотрим результаты моделирования для двух стадий выращивания монокристалла кремния диаметром 20 см и длиной 100 см, соответствующих 10 и 30% длины цилиндрической части слитка L_k. В этом случае предполагалось отсутствие теплового экрана 8, показанного на схеме теплового узла на рис. 1.

На рис. 2,a приведена ростовая стадия, соответствующая 10%-й длине выращенной цилиндрической части слитка. Скорости вращения кристалла и тигля соответствовали следующим параметрам: Re = 6 × 10⁴, Ros = –0.446, Gr = 2 × 10⁹. Структура течения является многовихревой и изменяющейся во времени. Дробление течения на ряд вихрей приводит к снижению теплопереноса от стенок тигля к ФК. Наличие в центре тигля направленного ко дну течения является причиной W-образного профиля ФК.

Рис. 2

На рис. 2,b показана картина теплопереноса в расплаве на стадии, соответствующей 30%-ной длине цилиндрической части слитка. По сравнению с предыдущим случаем скорость вращения тигля была несколько увеличена (до числа Ros = 0.54). Подъем тигля относительно нагревателя приводит к уменьшению влияния бокового нагрева и снижению числа Грасгофа до значения 1.7 × 10⁹. В этом случае происходит укрупнение вихря, прилегающего к боковой стенке тигля, за счет поглощения более мелких вихрей, что способствует лучшему транспорту тепла от тигля к ФК и более равномерному прогреву центральной области расплава. Такие изменения вызывают формирование вогнутой формы ФК.

Результаты расчета теплового поля в кристалле используются в последующем расчете процессов дефектообразования на основе построения “тепловой истории” для конкретного процесса выращивания кристалла кремния. Это предполагает получение аналитических выражений для теплового поля в кристалле T(t, r, z) в любой момент времени и для любой точки растущего кристалла с учетом формы ФК. С этой целью и вследствие того, что тепловые процессы достаточно быстрые, весь ростовой процесс представляется в виде последовательности N независимых стационарных стадий. Изменяющимися в ростовом процессе геометрическими параметрами, оказывающими существенное влияние на теплообмен в тепловом узле и на изменение теплового поля в кристалле, являются положение тигельной сборки, высота расплава в тигле и длина кристалла. Изменяющимися по заданным программам являются также мощность нагревателя и скорость вытягивания кристалла V.

Методический подход [5] позволяет запомнить и хранить температурные поля (“тепловую историю” роста монокристалла) для всего процесса выращивания. Сущность его состоит в представлении температурных распределений вдоль оси и кромки кристалла в виде полиномов высокой (6–8-й) степени. Они рассчитываются по температурным полям в кристалле для каждой расчетной ростовой стадии (всего для N ∼ 20 стадий). Далее в математических моделях для расчета дефектообразования эти полиномы задаются как формульные параметры (см. следующий раздел).

Распределение температуры в растущем монокристалле играет важную роль в процессах дефектообразования. Поэтому в стандартную конструкцию теплового узла для метода Чохральского добавляется специальный тепловой экран 8, показанный на рис. 1. С помощью этого экрана можно воздействовать на температурное поле в растущем монокристалле. Существенное воздействие теплового экрана на температурное поле в растущем монокристалле диаметром 20 см и длиной цилиндрической части 60 см можно видеть из сравнения картин изотерм в монокристалле при его выращивании без экрана (a) и с экраном (b), приведенных на рис. 3.

Рис. 3

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

1 – вблизи ФК при 1683 K происходит рекомбинация точечных дефектов (вакансий ${v}$ и междоузельных атомов i), температурные градиенты вблизи ФК определяют остаточную концентрацию вакансий и междоузельных атомов после процесса их рекомбинации;

2 – при 1373–1323 K образуются микропоры и кластеры междоузельных атомов кремния, формируется размер и концентрация ростовых микродефектов;

3 – при 1223–1023 K происходит распад твердого раствора кислорода с миграцией центров зарождения оксидных дефектов (преципитатов);

4 – при 1023–723 K растут размеры преципитатов и формируется окончательная картина микродефектов в монокристалле.

3. Моделирование процессов дефектообразования в монокристаллах и пластинах кремния. Следующим шагом сопряженного математического моделирования является использование “тепловой истории” выращиваемого монокристалла для расчета процессов дефектообразования в различных температурных зонах растущего монокристалла кремния (см. п. 2). Рассмотрим последовательность этих расчетов на примере “тепловой истории” выращивания монокристаллов кремния диаметром 150 мм.

В температурной зоне 1 происходит процесс рекомбинации собственных точечных дефектов (см. п. 2). Аннигиляция вакансий и междоузельных атомов вблизи ФК в значительной степени определяет их последующую миграцию и распределение в кристалле. Скорость аннигиляции равна разности между скоростями прямой реакции рекомбинации и обратной реакции термической генерации пар собственных точечных дефектов. Скорость прямой (бинарной) реакции рекомбинации пропорциональна произведению концентраций точечных дефектов при их взаимодействии. Скорость обратной реакции связана со скоростью прямой реакции законом действующих масс. В выращиваемом по методу Чохральского кристалле кремния перенос вакансий и междоузельных атомов осуществляется в основном конвективной диффузией, процесс термодиффузии пренебрегается [6].

В качестве примера на рис. 4 рассмотрен участок выращиваемого монокристалла вблизи фронта кристаллизации длиной l = 23 см и радиусом r = 7.5 см, для которого приведены изолинии ${{C}_{{i{v}}}}$ = C_i – ${{C}_{{v}}}$ [10^–12 см^–3] остаточных собственных точечных дефектов (где C_i, ${{C}_{{v}}}$ – концентрации междоузельных атомов и вакансий) после их рекомбинации в монокристалле кремния в двух режимах: a – вакансионном: 1) – 20.0, 2) – 16.0, 3) – 4.5, 4) – 0.1, 5) 2.7; b – смешанном (вакансионно-междоузельном): 1) – 16.0, 2) – 7.8, 3) – 4.5, 4) – 0.1, 5) 8.6. Различия вызваны разным соотношением скорости вытягивания V и осевого температурного градиента на ФК.

Рис. 4

В случае (a) концентрация остаточных вакансий – значительная, так как скорость вытягивания высокая. Однако при ее снижении растет концентрация междоузельных атомов и при некотором значении V происходит сужение вакансионной области к центру кристалла с расширением междоузельной области, примыкающей к его боковой поверхности (b). Радиальные распределения концентрации вакансий и междоузельных атомов на верхней границе расчетной области характеризуют их остаточные значения, которые являются начальными данными для расчета процесса их агломерации в микродефекты в более холодной температурной зоне 2.

В температурной зоне 2 для моделирования процесса образования ростовых микродефектов требуются как скорости охлаждения слитка, так и концентрации остаточных собственных точечных дефектов после рекомбинации. Математическая модель [5] учитывает двухмерность теплового поля в кристалле и диффузию вакансий.

На рис. 5 показаны радиальные распределения плотности ${{N}_{{v}}}$, N_p [см^–3] и размеров микропор ${{R}_{{v}}}$, R_p [нм] (a) и оксидных частиц (b) для растущего кристалла радиусом r = 7.5 см и длиной l = 0.3L_k, выращиваемого в смешанном вакансионно-междоузельном режиме (здесь сплошные линии соответствуют плотности, пунктир – размерам дефектов).

Рис. 5

Максимальная концентрация микропор достигает значения ${{N}_{{v}}}$ = 5 × 10⁴ см^–3, их размер ${{R}_{{v}}}$ = 23 нм в пределах центрального кольца радиусом 3 см (a). Максимальная концентрация оксидных частиц достигает большего значения: N_p = 10⁷ см^–3, а их размер существенно меньше: R_p = 4.2 нм в пределах центрального кольца радиусом 4 см. Моделирование образования междоузельных кластеров проводится по аналогичной схеме. Они расположены вблизи боковой поверхности слитка кремния. Оценки показывают, что значения их максимальной плотности и радиуса составляют примерно 8 × 10⁴ см^–3 и 54 нм.

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

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

Изменение скорости вытягивания монокристалла во время ростового процесса существенным образом влияет на распределение микродефектов. На рис. 6 показано изменение скорости вытягивания V [мм/мин] в виде зависимости от относительной длины выращенного кристалла l/Lк (a), которое согласуется с распределением расчетных изолиний микродефектов (b) и экспериментальной картой “времени жизни” электрических носителей (c). Это можно заметить по совпадению расчетного и измеренного контуров границы вакансионных и междоузельных микродефектов (см. пунктирные стрелки, указывающие на эту границу). Пространственные сужения и расширения контура этой границы согласованы с изменением скорости вытягивания.

Рис. 6

Верификация результатов математического моделирования также проводилась для пластин кремния, вырезаемых из выращенных монокристаллов [7]. Эти пластины проходят высокотемпературную термообработку в специальных термических печах с ламповым нагревом. Процесс состоит из быстрого разогрева пластины до 1523 K, кратковременную выдержку (10–15 с) пластины при этой температуре, а затем быстрое ее охлаждение до 823 K. В результате термообработки достигается снижение объемной неоднородности и уровня микродефектов. Применение рассмотренной выше сопряженной математической модели позволило рассчитать пространственно-временные изменения температурного поля в пластине кремния и в соответствии с этими данными рассчитать процесс диффузии и рекомбинации точечных дефектов. Для верификации результатов расчетов была измерена объемная плотность кислородных преципитатов, по данным которой построен профиль концентрации вакансий ${{C}_{{v}}}$ по толщине пластины. Его сравнение с аналогичным расчетным профилем концентрации вакансий показало хорошее соответствие [8].

4. Моделирование напряженного состояния в монокристаллах кремния. Наличие сверхбольших напряжений представляет актуальную проблему для технологии выращивания монокристаллов кремния диаметром 40 см и длиной более 100 см, для решения которой предлагаются специальные поддерживающие устройства в верхней конусной области, где находится крепление кристаллического слитка [9].

Однако при меньших диаметрах наиболее актуальным является анализ термонапряженного состояния выращиваемого монокристалла, в частности, в областях наибольших осевых и радиальных градиентов, возникающих вблизи ФК. Для монокристалла кремния диаметром 10 см и длиной цилиндрической части 40 см расчеты проведены в термоупругом приближении. На рис. 7 приведены картины изолиний напряжений по критерию Мизеса рассчитанные в соответствии с параметрами, приведенными в табл. 3. Их максимальные значения – 6.5 МПа достигаются на боковой поверхности кристалла вблизи ФК, что обусловлено наличием наибольших градиентов температуры в этой области. При удалении от ФК напряжения быстро уменьшаются.

Рис. 7

5. Моделирование напряженного состояния в пластинах кремния. Оптимизация режимов термообработки пластин кремния позволяет избежать появления дислокаций и дефектов структуры также в постростовой период. Например, результаты расчета трехмерного напряженного состояния пластин кремния становятся весьма актуальными при их большом диаметре ∼20 и 30 см.

Для процесса высокотемпературной обработки при горизонтальном креплении тонких пластин (толщиной ∼1000 мк) на трех опорах такие расчеты были выполнены с учетом кристаллографической ориентации (рис. 8,a: 1 – пластина кремния, 2 – крепление пластины на кольцевой опоре).

Рис. 8

В работе [10] было проанализировано влияние формы и материала опор пластины, а также их пространственного расположения на сдвиговые напряжения по возможным пространственным системам скольжения дислокаций в кристаллической решетке кремния. Для выявления дислокационно-опасных плоскостей скольжения рассчитывался прогиб пластины вдоль действия гравитационных сил при заданной расстановке и толщине опор.

Установлено, что напряжения по критерию Мизеса σ [МПа] достигают максимальных значений в области контакта пластины с опорой (рис. 8,b), причем пластины наименее напряжены при расположении опор на расстоянии 0.7R от центра пластины, где R – радиус пластины. Для пластин диаметром 20 см значения максимальных напряжений меньше примерно в 2 раза, чем для пластин большего диаметра 30 см.

Увеличение диаметра пластины до 30 см вызывает увеличение значений максимальных деформаций примерно в 4.5 раза. Для этого случая на рис. 9 приведены графики: максимальной величины напряжения σ [МПа] и величины изгиба пластины δ [мк] в зависимости от ширины кольцевой опоры h [мк] на расстоянии 0.7R от центра пластины.

Рис. 9

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

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

Работа выполнена на вычислительной базе ИПМех РАН (тема № АААА-А20-120011690136-2) при поддержке гранта РФФИ 18-02-00036.