Напряжения в пленках диоксида кремния, напыленных из диэлектрических мишеней: результаты атомистического моделирования Текст научной статьи по специальности «Физика»

- ФИЗИКА КОНДЕНСИРОВАННОГО СОСТОЯНИЯ ВЕЩЕСТВА -

Напряжения в пленках диоксида кремния, напыленных из диэлектрических мишеней: результаты атомистического моделирования

Ф.В. Григорьев,1, * В.Б. Сулимов,1 А.В. Тихонравов1

1 Московский государственный университет имени М. В. Ломоносова, научно-исследовательский вычислительный центр Россия, 119991, Москва, Ленинские горы, д. 1, стр. 4 (Поступила в редакцию 09.11.2023; после доработки 22.11.2023; подписана в печать 28.11.2023)

Предложенный ранее метод молекулярно-динамического моделирования напыления тонких пленок из металлических мишеней адаптирован к случаю диэлектрических мишеней и применен к пленкам диоксида кремния. Возможность вылета из мишеней не только атомов кремния, но и кластеров с атомами кислорода учтена путем добавления в поток осаждаемых атомов молекул O=Si=O. Получены атомистические кластеры пленок при высокоэнергетическом и низкоэнергетическом напылении при различном процентном содержании молекул в потоке осаждаемых атомов. Рассчитаны величины компонент тензора напряжений. При высокоэнергетическом напылении наблюдаются сжимающие напряжения, при низкоэнергетическом — растягивающие. Абсолютные величины диагональных компонент тензора напряжений растут с увеличением доли молекул в потоке осаждаемых атомов.

PACS: 81.15.Aa УДК: 539.231

Ключевые слова: суперкомпьютерное моделирование, молекулярная динамика, диоксид кремния, тонкие пленки.

DOI: 10.55959/MSU0579-9392.79.2410503

ВВЕДЕНИЕ

Физическое напыление является одним из наиболее распространенных методов изготовления оптических покрытий различного назначения [1]. Наиболее широко распространенный пленкообразующий материал с низким показателем преломления — диоксид кремния. При напылении пленок диоксида кремния в качестве мишени могут использоваться как кремниевые мишени [2], так и мишени, изготовленные из стеклообразного БЮ2 (диэлектрические мишени) [3]. В обоих случаях атомные кластеры, переведенные из твердой мишени в газовую фазу испарением или распылением пучком высокоэнергетических ионов благородных элементов, распространяются в вакуумной камере от мишени к подложке и осаждаются на ней, образуя слои пленки [1]. Из кремниевой мишени в газовую фазу переходят только отдельные атомы и кластеры кремния, из диэлектрической в составе кластеров в газовой фазе будут также и молекулы различного состава. Такое различие в составе потоков осаждаемых на подложку частиц повлечет за собой различие в свойствах пленок.

Экспериментальное исследование влияния кластеров на свойства пленок затруднено, так как состав потока осаждаемых частиц зависит от многих параметров — энергетического и углового распределения ионов благородных элементов, взаимодействующих с мишенью, плотности потока ионов,

* E-mail: fedor.grigoriev@gmail.com

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

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

1. МЕТОД МОДЕЛИРОВАНИЯ

Схема напыления тонких пленок представлена на рис. 1. При атомистическом моделировании процесса роста пленок выделяется микроскопический кластер с характерным размером 10-100 нм и числом атомов от нескольких тысяч до нескольких миллионов. Этот кластер используется для расчета параметров пленки.

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

Рис. 1. Схема напыления тонких пленок и область мо-лекулярно-динамического моделирования

Моделирование проведено развитым ранее методом, основанном на МД [4]. Явным образом моделируется только последняя стадия процесса напыления, когда атомы из мишени достигают поверхности подложки и ранее напыленных слоев пленки. Возможность выхода из мишени кластеров с атомами кислорода учтена добавлением в поток частиц трехатомных молекул 0=81=0, обладающих той же стехиометрией, что и материал мишени и растущей пленки. Доля молекул в потоке атомов определяется как отношение числа молекул 0=81=0 к общему числу атомов кремния в потоке частиц:

f = N(O=Si=O)/N.

(1)

Моделирование напыления организовано как пошаговая процедура. В начале каждого шага атомы кислорода и кремния, а также молекулы 0=81=0 располагаются случайным образом в верхней части области моделирования. Ориентация молекул выбирается случайным образом. Если в результате такого размещения центры атомов оказываются ближе, чем 0.2 нм (за исключением атомов кислорода и кремния в составе молекулы 0=81=0), координаты одного из двух атомов выбираются снова случайным образом, пока указанное условие не выполнится. Величина начальной скорости атомов кремния соответствует кинетической энергии 10 эВ при высокоэнергетическом напылении

и 0.1 эВ при низкоэнергетическом напылении. Начальные скорости атомов кислорода, входящих в состав молекул 0=81=0, равны по величине и направлению скорости атомов кремния. Начальные скорости атомов кислорода, не входящих в состав молекул, соответствуют кинетической энергии 0.1 эВ как для высокоэнергетического напыления, так и для низкоэнергетического напыления. Направление начальных скоростей определяется выбранным углом осаждения и одинаково для всех частиц, поскольку, как сказано выше, в эксперименте размер подложки много меньше расстояния от мишени до подложки.

После того, как заданы начальные положения и начальные скорости всех осаждаемых частиц, запускается моделирование в МУТ-ансамбле (постоянные число частиц, объем и температура, равная комнатной). Для поддержания постоянной температуры используется термостат Берендсена [5]. Длительность одного цикла составляет 10 пс. Как было показано ранее [6], такой длительности достаточно, чтобы атомы достигли поверхности пленки, отдали большую часть начальной кинетической энергии атомам пленки и подложки и образовали с ними химические связи. По окончании цикла вертикальный размер кластера увеличивается на 0.01 нм, чтобы учесть рост толщины пленки, и цикл повторяется. Общее число циклов моделирования достигало двух тысяч.

Потенциальная энергия взаимодействия между атомами рассчитывалась с использованием парного эмпирического силового поля БЕ81Ь [4, 7]:

U.

QiQj/тц + Ai./r j - Bij/r., (2)

где 4i(j) — заряд i(j)-ro атома, qo = —0.5qsi = —0.65e, Aij и Bij — параметры потенциала Леннарда-Джонса, rij — расстояние между атомами, Agio = 4.6 х 10-8 Дж-(нм)12/моль,

\12

моль,

Asísí = Aoo = 1.5 х 10-6 кДж-(н BS;o = 4.2 х 10-3 кДж-(нм)6/моль, Dsísí = BOO = 5 х 10-5 кДж-(нм)6/моль. Это же силовое поле используется и для атомов в составе молекул O=Si=O.

Молекулярно-динамическое моделирование проводилось с использованием программы GROMACS [8] на суперкомпьютерном комплексе МГУ имени М. В. Ломоносова [9]. Для визуального анализа структур использовалась программа VMD (Visual Molecular Dynamics) [10].

2. РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

2.1. Моделирование напыления пленок и расчет напряжений в растущей пленке

Горизонтальные размеры подложки при моделировании напыления взяты равными 10 нм и 30 нм, как и в предыдущих численных экспериментах [11]. Подложка получена из кристаллической структуры путем ее нагрева и охлаждения, детали

90 80 70 60

>

50 40 30 20

Д(80 = 0.1эВ

СXX

ГТ

иуу

0 Ю 20 30 40 50

/

-100

-120

Е(51)=\0эВ

Ь -140

-160

-180

Рис. 2. Зависимость главных компонент тензора напряжений а (МПа) от доли 1 (%) молекул 0=Я1=0 в потоке осаждаемых атомов

процедуры описаны в [12]. Скорость изменения температуры как при нагреве, так и при охлаждении составила 2 К/пс, что является типичной величиной при молекулярно-динамическом моделировании [13]. Конечная толщина подложки составила 5.5 нм.

Компоненты тензора напряжений равны компонентам тензора давления, взятым со знаком минус. Компоненты тензора давления, усредненные по траектории МД-моделирования атомистического кластера напыленной пленки, рассчитывались по следующим выражениям [8, 11, 14]:

Р

(К - Е),

(3)

где V — объем области моделирования, К и Е — тензоры кинетической энергии и вириала соответственно:

К

1

N

1 ^^ п 1

N

Е

г,3=1,г=3

% ® Гу. (4)

Здесь символ «(Е>» означает прямое произведение двух векторов.

Результаты расчета диагональных компонент тензора напряжений приведены на рис. 2, 4. Наблюдаются как сжимающие напряжения (а < 0), так и растягивающие (а > 0). Сжимающие напряжения характерны для высокоэнергетического напыления при энергии атомов кремния Е(Б1) = 10 эВ. Знак напряжений и их величины находятся в интервале экспериментальных значений [15, 16].

Как видно, с ростом доли O=Si=O молекул в потоке осаждаемых частиц абсолютные величины диагональных компонент тензора напряжений возрастают (рис. 2, правая сторона). Это можно объяснить следующим образом. При высокоэнергетическом напылении осаждаемые частицы, попадая в подложку и ранее напыленные слои пленки, «расталкивают» соседние атомы, что приводит к росту давления на стенки области моделирования.

Вероятно, это давление возрастает с увеличением массы и размера высокоэнергетических частиц, что и происходит с ростом доли молекул O=Si=O в потоке частиц.

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

При низкоэнергетическом напылении в структуре пленки появляются наноразмерные поры (рис. 3, левая сторона), которые отсутствуют при высокоэнергетическом напылении (рис. 3, правая сторона). Образование пор, в том числе крупных, — случайный процесс и поэтому может быть причиной некоторых флуктуаций при расчете компонент тензора напряжений.

Помимо диагональных компонент тензора напряжений, рассчитывались также и недиагональные, аХу и аух. Поскольку тензор напряжения симметричен, должно выполняться равенство аху = аух, что подтверждено в расчетах. Абсолютные величины недиагональных компонент не превышают десятых долей МПа для всех параметров моделирования. Таким образом, касательные напряжения, возникающие при напылении пленок диоксида кремния, незначительны.

Была рассчитана зависимость главных компонент тензора напряжений от толщины пленки. Результаты представлены на рис. 4.

Во всех случаях значения ахх и ауу растут по абсолютной величине при увеличении толщины пленки вплоть до Н ~ 25 нм, а затем выходят на константу. Отмеченные выше тенденции в зависимости ахх и ауу от энергии напыления и доли молекул в потоке осаждаемых атомов сохраняются при росте толщины пленки.

£(80=0.1 эВ, /-30%

£(80=10 эВ, /=30%

Рис. 3. Различия в структуре атомистических кластеров при высокоэнергетическом и низкоэнергетическом напылении. Нижняя часть кластера, окрашенная желтым, соответствует подложке

а £(80 = 0.1эВ 100

80

Ь 60 40 20

Е(&1) = 10 эВ -80

-100

/= 50%

10

15

25

30

Ь -140

-160

-180

-200

/ = эи /о

10

15

20 А

25

30

Рис. 4. Зависимость главных компонент тензора напряжений а (МПа) от толщины пленки Н (нм), Е (Я1) — энергия напыляемых атомов, / (%) — доля молекул 0=Я1=0 в потоке осаждаемых атомов

2.2. Моделирование термически индуцированных напряжений

Возникновение и рост напряжений в результате нагрева пленок лазерным излучением высокой интенсивности в настоящее время рассматривается в качестве одной из причин разрушения оптических покрытий [17]. В представленной работе мы провели моделирование нагрева пленки и возникающих при этом напряжений в зависимости от энергии

напыления и доли молекул O=Si=O в потоке осаждаемых атомов.

Нагрев кластера, напыленного до толщины Н = 20 нм, моделировался в МУХ-ансамбле. Температура равномерно увеличивалась с 300 К до 1500 К в течение 500 пс. Такой режим разогрева соответствует импульсам лазера наносекундой длительности. Нагрев моделирует передачу энергии лазерного излучения, поглощенную электронной подсистемой пленки, ядерной подсистеме по безыз-

лучательным каналам за счет электрон-фононного взаимодействия. Конечная геометрия являлась стартовой для последующего моделирования в течение 100 пс при постоянной температуре 1500 К в NVT ансамбле. Диагональные компоненты тензора напряжений усреднялись по последней траектории длительность 100 пс. То, что после окончания притока тепла температуру нагретой области пленки можно считать постоянной, следует из оценки характерного времени теплопередачи из уравнения теплопроводности: т = рЬ2/х, где с = 103 Дж/(кг-К), р = 2.2 х 103 кг/м3, X = 2.0 Вт/(м-К) — теплоемкость, плотность и коэффициент теплопроводности стеклообразного диоксида кремния, Ь ~ 0.1 мкм — характерная толщина поглощающего слоя. Получаем т ~ 10 нс ^ 100 пс. Таким образом, приближение постоянной температуры при расчете напряжений является обоснованным.

Таблица. Диагональные компоненты тензора напряжений ахх и Оуу (МПа) при нагреве пленки, Е (эВ) — энергия атомов кремния, / (%) — доля молекул 0=Я1=0 в потоке осаждаемых атомов

E 0.1 10

f 0 50 0 50

& XX -43 -12 -220 -230

ayy -40 -9 -217 -239

При низкоэнергетическом напылении оба диагональных компонента тензора напряжений изменили знак и стали сжимающими. Такое изменение вызвано стремлением вещества пленки увеличить свой объем при нагреве. Видно, что термически индуцированные напряжения в пленках, напыленных потоком частиц с высоким содержанием молекул O=Si=O, несколько ниже по абсолютной величине, чем в пленках, напыленных поток атомов. При высокоэнергетическом напылении знак диагональных

компонент тензора напряжений остался прежним, однако абсолютные величины ахх и ауу выросли, причем в случае / = 0 это увеличение существенно больше, чем в случае f = 50%. При этом разность в величинах ахх и ауу при переходе от / = 0 к f = 50% не очень значительна.

ЗАКЛЮЧЕНИЕ

В работе проведено молекулярно-динамическое моделирование процесса напыления тонких пленок диоксида кремния из диэлектрических мишеней. Возможность выхода из мишеней атомистических кластеров учитывалась добавлением в поток осаждаемых частиц трехатомных молекул O=Si=O. Исследовались случаи высокоэнергетического и низкоэнергетического напыления.

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

Также рассчитана зависимость напряжений от толщины пленки. Показано, что во всех случаях абсолютные величины диагональных компонентов тензора напряжений монотонно растут, пока толщина пленки не достигает приблизительно 25 нм, а затем выходят на константу.

Исследованы термически индуцированные напряжения, вызванные нагревом пленки из-за поглощения лазерного излучения. В пленках, полученных низкоэнергетическим напылением, при нагреве до 1500 К за 0.5 нс напряжения становятся сжимающими. Абсолютные величины главных компонентов тензора напряжений ниже в пленках, напыленных потоком частиц с высоким содержанием молекул.

Работа выполнена при поддержке Российского Научного Фонда, грант 23-11-00047.

[1] Piegari A. Optical Thin Films and Coatings. F.F., Ed.; Elsevier, 2018;

[2] Bernstein J., Gerlach J.W., Finzel A., Bundesmann C. I. // Eur. Phys. J. B 95, 39 (2022).

[3] Telesh E. V., Dostanko A.P., Gurevich O. V. // J. Appl. Spectrosc. 85, 67 (2018).

[4] Григорьев Ф.В. // Вестн. Моск. ун-та. Физ. Астрон. 6, 93. (2015). (F. V. Grigoriev // Moscow Univ. Phys. Bull. 70, N 6. 521 (2015).)

[5] Berendsen H.J.C., Postma, J.P.M., Van Gunsteren W.F. et al. // J. Chem. Phys. 81, 3684 (1984).

[6] Григорьев Ф. В., Сулимов В. Б., Кондакова О. А. и др. Вестн. Моск. ун-та. Физ. Астрон. 3, 80. (2013). (F. V. Grogoriev, V. B. Sulimov, O. A. Kondakova et al. // Moscow Univ. Phys. Bull. 68, N 3, 259 (2013).)

[7] Grigoriev F. V., Sulimov A. V., Kochikov I. V. et al. .. In Proceedings of the Proceedings of SPIE -

The International Society for Optical Engineering; Bellingham, WA, United States, 9627, 962701 (2015).

[8] Abraham M.J., Murtola T., Schulz R. et al. GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1-2, 19-25

[9] Voevodin V. V., Antono, A.S., Nikitenko D.A., et al. // Supercomput. Front. Innov. 6, N 2 (2019).

[10] Humphrey W., Dalke A., Schulten K. // J. Mol. Graph. 14, 33 (1996).

[11] Grigoriev F. V., Sulimov V.B., Tikhonravov A. V. // Coatings 10, 220 (2020).

[12] Grigoriev F. V., Sulimov A. V., Katkov E. V. et al. // Appl. Opt. 56, C87 (2017).

[13] Kirova E.M., Norman G.E., Pisarev V. V. Russ. J. Phys. Chem. A 92, 1865 (2018).

[14] Al M., Webb III. E. // Nanomater. Nanotechnol. 8, 1

(2018).

[15] Begou T., Lumeau J. // Opt. Lett. 42, 3217 (2017).

[16] Abadias G., Chason E., Keckes J. et al. // J. Vac. Sci. Technol. A 36, 20801 (2018).

[17] Sherman B., Liou H.-C., Balogun O. //J. Appl. Phys. 118, 135303 (2015).

[18] Леко В.К., Мазурин О.В. Свойства кварцевого стекла. Ленинград: Наука, 1985.

Stresses in silicon dioxide films deposited from dielectric targets: results

of atomistic modeling

F.V. Grigoriev", V.B. Sulimov, A.V. Tikhonravov

Research Computing Center, Lomonosov Moscow State University. Moscow 119991, Russia

E-mail: afedor.grigoriev@gmail.com

The previously proposed method of molecular dynamics simulation of the deposition of thin films from metal targets is adapted to the case of dielectric targets and applied to silicon dioxide films. The possibility of not only silicon atoms leaving the target, but also clusters with oxygen atoms was taken into account by adding O=Si=O molecules to the flow of deposited atoms. Atomistic film clusters were obtained by high-energy and low-energy deposition at different percentages of molecules in the flow of deposited atoms. The values of the stress tensor components are calculated. With high-energy deposition, compressive stresses are observed, with low-energy deposition, tensile stresses are observed. The absolute values of the diagonal components of the stress tensor increase with increasing proportion of molecules in the flow of deposited atoms.

PACS: 81.15.Aa

Keywords: high-performance simulation, molecular dynamics, thin film growth, deposition process, silicon dioxide, porosity. Received 9 November 2023.

English version: Moscow University Physics Bulletin. 2024. 79, No. 1. Pp. 52-57. Сведения об авторах

1. Григорьев Федор Васильевич — доктор физ.-мат. наук, вед. науч. сотр.; тел.: (495) 939-32-53, e-mail: fedor.grigoriev@gmail.com.

2. Сулимов Владимир Борисович — доктор физ.-мат. наук, зав. лаб.; тел.: (495) 939-32-53, e-mail: vs@dimonta.com.

3. Тихонравов Александр Владимирович — доктор физ.-мат. наук, науч. руководитель НИВЦ МГУ им. М. В. Ломоносова; тел.: (495) 939-54-24, e-mail: tikh@srcc.msu.ru.