Термическое оксидирование – это процесс создания очень тонкой в бытовом понимании (десятки/сотни нанометров) пленки диоксида кремния $\rm Si O_2$ на поверхности кремния $\rm Si$ при высоких температурах. Данная технология является важнейшей частью процесса изготовления интегральных микросхем и МОП-транзисторов (металл-оксид-полупроводник), которые составляют основу современной электроники. Если поместить кремний $\rm Si$ в аэробную среду, то начнется процесс окисления $\rm Si + O_2 \rightarrow SiO_2$. Для описания кристаллического кремния и плоских пленок оксида достаточной толщины (начиная с десяти нанометров) хорошо себя показывает модель Дила-Гроува (1965), но детальное описание механизма роста в других условиях до сих пор остается нерешенной научной задачей.

Кремний является непрозрачным для видимого света, в то время как его диоксид прозрачен. Эти оптические свойства в купе с тем, что пленка диоксида имеет толщину десятки нанометров, создают идеальные условия для очень точного оптического контроля толщины слоя $\rm SiO_2$.

В рамках нашей задачи рост пленки будет исследоваться с помощью эллипсометрии - мощной методики исследования оптических свойств тонких пленок. Она основана на сравнении амплитудных коэффициентов отражения волн $s$-поляризации и $p$-поляризации (от нем. "senkrecht" – перпендикулярный, "parallel" – параллельный) от поверхности пленки. Результаты измерений записываются с помощью эллипсометрических углов $\Psi$, $\Delta$:

\[ \frac{r_p}{r_s} = \tan \Psi \cdot e^{i \Delta}, \]

где значение $\Psi$ ответственно за отношение амплитуд при отражении $p$- и $s$-поляризованных волн, а значение $\Delta$ – за разность фаз между ними.

При эллипсометрии измеряется именно относительное различие между волнами (отношение амплитуд и разность фаз) с разной поляризаций, поэтому она очень устойчива к различным флуктуациям в мощности источника излучения и чувствительности приемника. Также эллипсометрия позволяет измерять характеристики образца без непосредственного механического контакта с ним, что позволяет проводить измерения в реальном времени и in situ (с лат. – “на месте”, непосредственно).

В течение всей задачи длина волны падающего света обозначается $\lambda$. В оптическом диапазоне у веществ нет значительных магнитных свойств, поэтому у всех рассматриваемых материалов $\mu=1$. 

Внимание! Помещайте файлы, которые хотите приложить к отчету, в папку report.

Внимание! Для решения части B в части А достаточно сделать только пункт А7.

Часть А. Теоретические принципы эллипсометрии

Уравнения Максвелла в немагнитном диэлектрике имеют вид
\[
\begin{split}
\DeclareMathOperator{\Div}{div}
\DeclareMathOperator{\Rot}{rot}
\DeclareMathOperator{\Grad}{grad}
\Div {\vec{D}} = 0 \quad \quad & \Rot{\vec{E}} = - \frac{\partial \vec{B}}{\partial t} \\
\Div{\vec{B}} = 0 \quad \quad & \Rot{\vec{B}} = \mu_0 \frac{\partial \vec{D}}{\partial t}
\end{split}
\]
Для диэлектрика:
\[\vec{D} = \varepsilon_0 \varepsilon \vec{E}\]Из этих уравнений следует, что нормальные компоненты векторов $\vec{D}$ и $\vec{B}$ и направленные вдоль поверхности компоненты векторов $\vec{E}$ и $\vec{B}$ не изменяются на границе диэлектрик-диэлектрик.

В плоской волне волновой вектор, вектор напряженности электрического поля и вектор индукции магнитного поля $\vec{k}$, $\vec{E}$ и $\vec{B}$ взаимно перпендикулярны и образуют правую тройку. Также справедливы соотношения:
\[\vec{B} = \frac{1}{\omega} [\vec{k} \times \vec{E}],\]\[|\vec{k}| = \sqrt{\varepsilon} \omega/c,\]где $\omega$ – частота волны. Cкорость света $c = 1/\sqrt{\varepsilon_0 \mu_0}$.
В рамках этой задачи мы будем рассматривать взаимодействие плоской волны со слоистыми структурами, поэтому выработаем общий подход для теоретического описания такой системы. Пусть пленка толщиной $d$ из материала с диэлектрической проницаемостью $\varepsilon_2$ находится между двух полубесконечных диэлектриков с проницаемостями $\varepsilon_1$ и $\varepsilon_3$.

При падении плоской волны на пленку возникает отраженная волна в среде $\varepsilon_1$, прошедшая волна в среде $\varepsilon_3$ и суперпозиция волн идущих в разных направлениях в среде $\varepsilon_2$.
В $s$-поляризации направление электрического поля фиксировано, поэтому можно перейти к рассмотрению скалярных комплексных амплитуд:
\[
\begin{cases}
\tilde{E}_1 = \mathcal{E} e^{ik_{1x}x+ik_{1y}y} + r \cdot\mathcal{E} e^{-ik_{1x}x+ik_{1y}y} \\
\tilde{E}_2 = a \cdot \mathcal{E} e^{ik_{2x}x+ik_{2y}y} + b \cdot \mathcal{E} e^{-ik_{2x}x+ik_{2y}y} \\
\tilde{E}_3 = t \cdot \mathcal{E} e^{ik_{3x}(x-d)+ik_{3y}y}
\end{cases}
\]
В $p$-поляризации направление магнитного поля фиксировано, поэтому проще перейти к рассмотрению скалярных комплексных амплитуд уже магнитного поля:
\[
\begin{cases}
\tilde{B}_1 = \mathcal{B} e^{ik_{1x}x+ik_{1y}y} + r \cdot\mathcal{B} e^{-ik_{1x}x+ik_{1y}y} \\
\tilde{B}_2 = a \cdot \mathcal{B} e^{ik_{2x}x+ik_{2y}y} + b \cdot \mathcal{B} e^{-ik_{2x}x+ik_{2y}y} \\
\tilde{B}_3 = t \cdot \mathcal{B} e^{ik_{3x}(x-d)+ik_{3y}y}
\end{cases}
\]

A1  ^0.20 Покажите, что тангенциальные компоненты волновых векторов не изменяются на границе диэлектрик-диэлектрик: $k_{1y}=k_{2y}=k_{3y}$.

А2  ^0.40 Найдите выражения для $k_{1x}$, $k_{2x}$ и $k_{3x}$. Выразите через $\omega$, $c$, $\theta_1$, $\varepsilon_1$, $\varepsilon_2$, $\varepsilon_3$.

A3  ^0.60 Для волны с $s$ поляризацией (электрическое поле направлено вдоль поверхности) из граничных условий на тангенциальные компоненты полей получите систему линейных уравнений на коэффициенты $r$, $a$, $b$ и $t$.

A4  ^0.60 Найдите коэффициент отражения. Приведите ответ к виду:
\[ r_s = \frac{A +B e^{2i k_{2x}d}}{AB + e^{2i k_{2x}d}}.\]Вещественные числа $A$, $B$ выразите через $k_{1x}$, $k_{2x}$, $k_{3x}$.

Для удобства введем переменные $\kappa_{ix}=\cfrac{k_{ix}}{\varepsilon_i}$ для $i\in[1,3].$

A5  ^0.60 Для волны с $p$ поляризацией (магнитное поле направлено вдоль поверхности) из граничных условий на тангенциальные компоненты полей получите систему линейных уравнений на коэффициенты $r$, $a$, $b$ и $t$.

A6  ^0.60 Найдите коэффициент отражения. Приведите ответ к виду:
\[ r_p = \frac{A +B e^{2i k_{2x}d}}{AB + e^{2i k_{2x}d}}.\]Вещественные числа $A$, $B$ выразите через $\kappa_{1x}$, $\kappa_{2x}$, $\kappa_{3x}$.

Проведем предварительные расчёты примерного вида эллипсометрического спектра тонкого слоя диоксида кремния на подложке из кремния. До конца части А будем пренебрегать оптической дисперсией, и будем считать $\varepsilon_\text{Air}=1.00$, $\varepsilon_{\rm Si} = 17.2 - 0.430 i$, $\varepsilon_{\rm SiO_2}=2.13$. Заметим, что кремний обладает комплексной составляющей диэлектрической проницаемости, ответственной за поглощение, и поэтому он является непрозрачным. Использование комплексных диэлектрических проницаемостей никак не меняет вид полученных выражений, если определить тригонометрические функции комплексного параметра.

C помощью программы A.py вы можете строить графики $\Psi$ и $\Delta$ для разных углов падения $\theta_1$ и толщин $d$. На вход программа получает название серии Name, угол падения $\theta_1$ из воздуха в градусах и толщину $d$ пленки диоксида виде таблицы, записанной в файле A_in.txt.

Пример заполнения таблицы в файле A_in.txt:

Для входных данных заложены ограничения $20\,^\circ \leq \theta_1 \leq 80\,^\circ$, $d \leq 250~\text{нм}$. Название эксперимента не должно содержать пробелов и табуляции. Разделитель между целой и дробной частью чисел – точка. Между столбцами разделитель –табуляция.

A7  ^1.50 С помощью программы A.py постройте графики зависимости $\Psi$ и $\Delta$ от $\lambda$ при разных углах падения $\theta_1$ и толщинах $d$. Проведите количественное исследование и определите, при каком угле падения $\theta_\text{opt}$ эллипсометрический угол $\Psi$ наиболее чувствителен (наибольшее значение $|\Delta \Psi/\Delta d|$) к изменениям $d$ в диапазоне от $10 нм$ до $50нм$. Найдите значение $\theta_\text{opt}$ с точностью до градуса.

Начиная с пункта B3, для повышения точности используйте $\theta_{opt}$ в качестве $\Phi$ .

Часть В. Кинетика оксидирования

В данной части угол падения луча на пленку будет обозначаться $\Phi$.

Вы вместе с коллегой занимаетесь изучением образования тонких оксидных пленок на поверхности кремния при его длительном прокаливании на сухом воздухе. Реальный кремний $\rm Si$ и диоксид кремния $\rm SiO_2$ обладают оптической дисперсией, поэтому графики $\Psi, \Delta$ от $\lambda$, которые будут предложены вам в этой части, могут отличаться от тех, которые вы строили в рамках пункта А7.

Внимание! Оптическая дисперсия $\rm Si$ и $\rm Si O_2$ автоматически заложена во всех программах, с которыми вы будете работать в этой части.

Ваш коллега вчера допоздна работал на эллипсометре и забыл записать значение угла $\Phi$. Он измерял толщину оксидной пленки на двух образцах и результаты измерений передал вам в виде файлов B1.txt и B2.txt.

B1  ^0.70 С помощью программы B1.py подберите вручную значение $\Phi$ и толщины $d$ пленки $\rm SiO_2$. Программа B1.py автоматически считывает значения из файла B1.txt и выводит аппроксимирующий график $\Psi,\Delta$ от $\lambda$, для введенных вами значений угла $\Phi$ в градусах и толщины $d$ пленки диоксида в нанометрах. Для входных данных заложены ограничения $20\,^\circ \leq \Phi \leq 80\,^\circ$, $d \leq 250~\text{нм}$.

B2  ^0.70 С помощью программы B2.py подберите вручную значение $\Phi$ и толщины $d$ пленки $\rm SiO_2$. Программа работает аналогично предыдущему пункту.

Для подробного изучения скорости окисления вы можете заказывать у вашего коллеги эллипсометрические данные, измеренные с окисляющейся в реальном времени подложки кремния. В начальный момент времени $t=0$ он помещает чистый кремний внутрь печки с температурой $T$ и раз в $2~\text{мин}$ снимает спектр $\Psi, \Delta$ от $\lambda$ при фиксированном угле $\Phi$, который вы ему сообщите.

В файле B_in.txt вы задаете название Name эксперимента, угол $\Phi$ в градусах, $T$ в градусах Цельсия и длительность $t$ эксперимента в минутах в виде таблицы:

Для входных данных заложены ограничения $20\,^\circ \leq \Phi \leq 80\,^\circ$, $800\,^\circ\text{C} \leq T \leq 1200\,^\circ\text{C}$, $t \leq 200~\text{мин}$. Название эксперимента не должно содержать пробелов и табуляции. Разделитель между целой и дробной частью чисел – точка. Между столбцами разделитель – табуляция. Далее вы запускаете программу gen_data.exe которая в ускоренном режиме симулирует проведение эксперимента вашим коллегой. Вы можете прервать выполнение этой программы в любой момент времени и результаты уже выполненных измерений сохранятся у вас на компьютере. Вы можете автоматически обработать данные измерений с помощью программы B.py, которая использует тот же файл B_in.txt. Эта программа на основании результатов работы gen_data.py подбирает значение $d$ для каждой точки и рассчитывает статистическую погрешность этого значения. Имейте в виду, что обработка может занять несколько минут. Результаты работы программы B.py находятся в папке /res в виде таблиц и графиков $d(t)$. Пример автоматически сгенерированной программой B.py таблицы в папке /res:

B3  ^0.30 Получите данные для зависимости толщины пленки $d$ диоксида кремния от времени $t$ при окислении при температуре $t=1000\,^\circ\text{C}$ в течение $120~\text{мин}$. Приложите полученный с помощью программы B.py график к отчету.

Теперь рассмотрим теоретические аспекты термической оксидации. Согласно модели Дила-Гроува за скорость увеличения толщины слоя оксида кремния ответственны три процесса:

1. диффузия кислорода к поверхности $\rm SiO_2$;
2. диффузия кислорода сквозь слой $\rm SiO_2$;
3. реакция окисления $\rm Si + O_2 \rightarrow SiO_2$ на границе раздела $\rm SiO_2$/$\rm Si$.

Таким образом, концентрация $C$ кислорода (размерность $\text{моль}/\text{м}^3$) изменяется при движении вглубь образца. В воздухе, достаточно далеко от подложки, концентрация кислорода равна $C^*$; внутри слоя диоксида кремния концентрация линейно спадает от $C_0$ до $C_i$; внутри кремния концентрация кислорода равна нулю.

Рассмотрим модель для описания указанных процессов. Удельные (т. е. на единицу площади поверхности) скорости протекания каждого из трех процессов $F_1$, $F_2$, $F_3$ (единицы измерения – $\text{моль}/(\text{с}\cdot\text{м}^2)$) будем определять следующим образом:

1. скорость диффузии $F_1 = h(C^* - C_0)$, где $h$ – некоторый коэффициент, характеризующий диффузию в газовой фазе.
2. Считая скорость диффузии $F_2 = \frac{D}{d}(C_0 - C_i)$, где $D$ – коэффициент диффузии кислорода в диоксиде, зависящий от температуры $T$,  и $d$ – толщина пленки диоксида кремния.
3. Скорость реакции окисления $F_3 = k C_i$, где $k$ – константа скорости реакции окисления, зависящая от температуры $T$.
    

B4  ^0.30 Выразите скорость $v$ изменения толщины пленки через скорость реакции окисления $F_3$, молярную массу $\mu_{\rm SiO_2}$ и плотность $\rho_{\rm SiO_2}$ диоксида кремния. Считайте, что граница воздух-диоксид неподвижна.

B5  ^0.50 Считая задачу псевдостационарной (времена установления всех процессов гораздо меньше времени нарастания пленки), выразите $C_0$ и $C_i$ через $C^*$, константы скорости $h$, $k$, $D$ и толщину пленки $d$.

B6  ^0.50 Получите дифференциальное уравнение для толщины $d$ пленки диоксида и получите его решение для начальных условий $d(t_0)=d_0$ в виде
\[ d^2 + Ad = B(t + \tau),\]где $\tau$ содержит в себе информацию о начальных условиях, а константы $A$ и $B$ связаны с кинетикой окисления. Выразите $A$ и $B$ через $h$, $k$, $D$, $C^*$, $\mu_{\rm SiO_2}$ и $\rho_{\rm SiO_2}$.

Экспериментальные исследования показывают, что при продувании воздухом (то есть при изменении профиля $C(x)$ над поверхностью) скорость процесса не изменяется, поэтому будем считать $h \gg k$.

В начале процесса окисления стационарная модель не применима, поэтому ее применяют начиная с минимальной толщины $d_0=25~\text{нм}$. Введите в программу B7.py данные из файла B_in.txt и определите параметры $B$ и $B/A$ автоматическим образом. Если максимальная толщина пленки, полученная в эксперименте, меньше $d_0$, то программа B7.py выведет "Final layer is too thin to approximate".

В случае если программа выдает ошибку "FileNotFoundError", последовательно запустите файлы gen_data.exe и B.py с нужными вам входными данными, предварительно занесёнными в B_in.txt.

B7  ^1.00 Получите значения параметров $B$ и $B/A$ для пяти разных температур $T$. Результаты занесите в файл B8_in.txt в виде таблицы.

Зависимость константы скорости реакции $k$ от температуры $T$ связывают с энергией $\Delta E_b$ разрыва связи $\rm Si - Si$ в кристалле кремния, и тогда связь между ними:
\[ k \propto e^{-\Delta E_b/k_BT}.\]Коэффициент диффузии в твердых веществах также связан с температурой через т.н. энергию активации $\Delta E_a$, и количественная зависимость имеет вид:
\[ D \propto e^{-\Delta E_a/k_BT} \]

С помощью программы B8.py вы можете построить графики $\ln(B \cdot \text{мин}/\text{нм}^2)$ от $1/T$ и $\ln(B/A \cdot \text{мин}/\text{нм})$ от $1/T$ и определить параметры автоматически проведенных прямых.

Постоянная Больцмана $k_B=1.38 \cdot 10^{-23}~\text{Дж}/\text{К}$, элементарный заряд $e=1.60\cdot 10^{-19}~\text{Кл}$

B8  ^1.50 Определите значения $\Delta E_b$ и $\Delta E_a$ в электрон-вольтах. Приложите линеаризованные графики к отчету.