CC BY
8
УДК 538.913+534-6 DOI: 10.15350/17270529.2020.3.27
АТОМАРНЫЙ МЕХАНИЗМ РАЗРУШЕНИЯ ТОНКОЙ ПЛЕНКИ ПРИ ВЫСОКИХ ТЕМПЕРАТУРАХ. МЕТОД МОЛЕКУЛЯРНОЙ ДИНАМИКИ
САЛАМАТОВ Е. И., ДОЛГУШЕВА Е. Б.
Удмуртский федеральный исследовательский центр Уральского отделения РАН, 426067, г. Ижевск, ул. Т. Барамзиной, 34
АННОТАЦИЯ. Методом молекулярной динамики на атомарном уровне исследован процесс разрушения ультратонкой (5 атомных слоев) алюминиевой пленки с ГЦК структурой. Расчеты проводились с использованием потенциала межатомного взаимодействия, аттестованного путем сравнения численных результатов с аналитическими, полученными в рамках теории упругости. Показано, что при температурах меньших 600 K система остается стабильной на временном интервале 0.6 ш, что позволяет анализировать низкочастотную часть спектра до частоты Утт = 0.0166 THz. Установлено, что определенная при T = 50 K частота колебаний самой длинной, для данной геометрии задачи, изгибной волны у0 = 0.033 THz, уменьшается с ростом температуры, соответственно, её период растёт. Как только период колебаний этой волны, названной нами «заторможенной модой», становится сравним с временем моделирования, её амплитуда начинает расти. Показано, что процесс разрушения пленки стартует после достижения амплитудой «заторможенной моды» некоторого критического значения.
КЛЮЧЕВЫЕ СЛОВА: метод молекулярной динамики, преобразование Фурье, закон дисперсии колебаний тонкой пленки, «заторможенная мода».
ВВЕДЕНИЕ
Применение в промышленности метастабильных нанопленок, обладающими уникальными свойствами, в значительной мере затруднено тем обстоятельством, что в свободном состоянии они нестабильны [1 - 4]. В связи с этим, первоочередной задачей является повышение стабильности пленок. Решение этой проблемы требует детального понимания условий и факторов, которые играют решающую роль в процессе разрушения пленки, определения ключевых изменений в системе, однозначно указывающих на начало этого процесса.
Из-за очень малого времени жизни свободных кристаллических пленок проведение экспериментальных исследований в настоящее время довольно проблематично. Для решения задач, связанных с рассмотрением на атомарном уровне динамических процессов, происходящих за малые временные интервалы в системах с большим количеством частиц, наиболее перспективным является метод молекулярной динамики (ММД).
Хорошо известно, что результаты, полученные с помощью ММД, зависят от используемого потенциала межатомного взаимодействия, геометрии и размера расчетной ячейки, выбора граничных условий, а также времени счета. Поэтому необходимо провести предварительную аттестацию потенциала, сравнивая рассчитанные свойства с известными физическими характеристиками моделируемых объектов, чтобы убедиться в достоверности описания. Геометрия расчетной ячейки должна позволять решить поставленную задачу, кроме того, её размер определяет максимальную длину волны Х™ах и, следовательно, минимальную частоту утщ, реализующиеся в системе. Это позволяет установить необходимое время счета ^ для того, чтобы зафиксировать пик в колебательном спектре, соответствующий минимальной частоте, и исследовать температурную зависимость утш.
Ранее нами была предложена методика тестирования результатов ММД, которая показала, что выбранный потенциал позволяет адекватно описать динамику сверхтонких пленок алюминия [5].
Таким образом, используя наши предыдущие результаты, в настоящей работе определены причины развития динамической неустойчивости тонкой пленки с ростом температуры, приводящей к ее разрушению, что позволило предложить разумный сценарий этого процесса на атомном уровне. Заметим, что при периодических граничных условиях в двух направлениях, необходимых для расчета динамических характеристик, невозможно добиться реального разрушения пленки: разрыва или скручивания. Поэтому под разрушением мы будем понимать такие изменения в системе, которые приводят к остановке программы, предполагая, что эти же изменения происходят в реальном случае на начальной стадии процесса разрушения.
МЕТОДЫ ВЫЧИСЛЕНИЙ (ИСХОДНЫЕ ДАННЫЕ)
Как уже отмечалось во введении, в настоящей работе мы будем отталкиваться от результатов, полученных в [5]. В ней исследовалась стабильность пленки ГЦК-А1 толщиной в 5 атомных слоев, базовая расчетная ячейка которой изображена на рис. 1. Вдоль осей у и г задавались периодические граничные условия, вдоль оси х - свободные.
Рис. 1. Модельный базовый кристаллит, имитирующий пленку.
Толщина пленки N = 5 моноатомных слоев Л1, N = 30, N2 = 84. Всего 6300 атомов
Проблема тестирования результатов динамических свойств, определенных с помощью ММД для тонких пленок, была решена путем сравнения вычисленного закона дисперсии с аналитическим, полученным в рамках теории упругости. Численные расчеты опирались на предварительно проведенный симметрийный анализ колебаний атомов пленки, что позволяет независимо вычислять частотные спектры для симметризованных координат, соответствующих изгибным ^х,ь) и сжимающим колебаниям [6] (рис. 2). Такой анализ является очень полезным при описании физических свойств пленок, поскольку колебания разной симметрии ответственны за различные свойства системы [7, 8].
В работе [5] найдены температуры (< 600 К), при которых система остается стабильной на временном интервале 0.6 м. Это позволяет анализировать низкочастотную часть спектра до минимальной частоты уШт = 0.0166 ТН и исследовать температурную зависимость колебаний самой длинной, для данной геометрии задачи, изгибной волны с частотой у0 = 0.033 ТНг, полученной при Т = 50 К.
О 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
V, ТНг
Рис. 2. Низкочастотный спектр изгибных ^х,ь) и сжимающих ^х,с) колебаний поверхностных атомов пленки, построенный при Т = 50 К.
Спектр gx,c увеличен в 30 раз. Вертикальной штриховой линией показано положение минимальной частоты = 0.0166 ТИг
РАЗРУШЕНИЕ ПЛЕНКИ
При приближении температуры к 800 ^ время жизни стабильного состояния пленки становится меньше 0.6 м. Это проявляется в быстром, экспоненциальном уменьшении размера базового кристаллита Ь2{() вдоль направления г при временах больших некоторого критического значения Заметим, что время ^ является функцией не только температуры, оно также зависит от способа выведения системы на данную температуру. Зависимости размера кристаллита Ь2{() от времени, приведенные на рис. 3, стартовали после выдержки системы в течение 0.1 м при 300 К, что позволяет считать условия достижения высоких температур (>770 ^ примерно одинаковыми для всех кривых Ь2(1).
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1: , Г18
Рис. 3. Временная зависимость размера кристаллита вдоль оси г, полученная при различных температурах. Штриховой линией показана Ь^) при Т = 785 К
Кроме размеров кристаллита отслеживалось изменение средних значений элементов тензора внутренних напряжений [9]:
(a)aß=~ ^
Л
V 2 i(zV j*'
(1)
Здесь индексы а, ß принимают значения х, y, z; V - объем кристаллита; i, j - номера атомов; Fj - силы и Ту - расстояние между атомами i и j.
При расчете элементов тензора учитывалось только слагаемое, возникающее непосредственно из-за сил, вириальный вклад не учитывался. Временные зависимости элементов тензора внутренних напряжений охх, ozz системы тоже претерпевают качественное изменение при временах близких к t0(T). Все полученные кривые хорошо аппроксимируются одним выражением:
Ф(0 = f + at + b• exp((t-to)/x). (2)
Параметры f0, a, b, t0, т как для размера кристаллита (L), так и для элементов тензора внутренних напряжений (о) при различных температурах были найдены путем подгонки с использованием нелинейного алгоритма наименьших квадратов (NLS) [10].
Так как временные характеристики, рассчитанные при температурах T > 770 K, подобны друг другу, остановимся более подробно на результатах, полученных для температуры 785 К. На рис. 4 показаны временные зависимости, полученные при T = 785 K для размера пленки вдоль оси z Lz(t) - (а) и диагональных элементов тензора напряжений Охх и Ozz - (b). Здесь же приведены их аппроксимации по формуле (2). На рисунках указано значение только одного подгоночного параметра - t0, определяющего момент времени качественного изменения соответствующих характеристик: для Lz t0 = 0.321 ns, а для о - t0 = 0.565 ns.
■о> о
ь
170 160 150 140 130 120 110
t0=0.321ns
а)
О 0.1 0.2 0.3 0.4 0.5
t, ns
^zz
Рис. 4. Временная зависимость: а) - размера базового кристаллита вдоль оси г Ь) - диагональных элементов тензора напряжения охх агг при Т = 785 К
Из сравнения этих графиков видно, что при заданной температуре первые признаки неустойчивости системы появляются раньше на кривой Ь2(().
На рис. 5 приведены мгновенные проекции всех атомов пленки на плоскость Х2 при Т = 785 К в различные моменты времени, значения которых указаны справа для каждого состояния. По этому рисунку можно проследить, как с течением времени изменяется форма пленки при приближении к нестабильности. Отметим, что мгновенная форма пленки совпадает с геометрией самой длинной низкочастотной изгибной волны (см. рис. 3, Ь из работы [5]), амплитуда которой растет со временем, что и обеспечивает линейное уменьшение Ь2 при временах I < ¿0. При временах больше ¿0 процесс ускоряется, происходит быстрое нарастание амплитуды изгибной волны, пленка изгибается, при этом длина волны резко уменьшается, соответственно, уменьшается Ь2 и, как только Ь2 становится равным трем межатомным расстояниям, счет прекращается.
Рис. 5. Мгновенные проекции всех атомов пленки на плоскость XI в различные моменты времени при Т = 785 К. Ьг(^) - изменение размера базового кристаллита вдоль оси г при этой же температуре
Проанализируем, как изменяется самая низкочастотная часть спектра системы с ростом температуры до 600 К. На рис. 6 представлены результаты расчета низкочастотного спектра изгибных колебаний, полученные при температурах: 50, 200, 300 и 600 К. Наблюдаемое с ростом температуры смягчение спектра вполне объяснимо - оно связано с уменьшением скорости поперечной звуковой волны в объемном образце (см. выражение (1) из [5]). При Т = 600 К период колебаний этой волны сравним с временем счета (0.6 м) и, очевидно, продолжает увеличиваться с ростом температуры. Частоты меньше величины не определяются из Фурье анализа, а колебательную моду, период которой больше времени счета, можно называть «заторможенной». Действительно, с самого начала расчета амплитуда самой длинной волны непрерывно растет, при этом длина (Хщах) её уменьшается, все это приводит к увеличению кривизны пленки (см. рис. 5). Если в момент времени ¿0 < ^ длина волны Хшах (а вместе с ней и длина кристаллита Ь2) достигает некоторого критического значения ((Ь2(0) - Ь2(^)/Ь2(0) ~ 0.01), то становится возможным флуктуационное притягивающее взаимодействие между атомами, расположенными на противоположных склонах гребня волны, на что указывает появление отрицательных значений элементов тензора о22 (см. рис. 4, Ь). После чего этот процесс развивается экспоненциально, происходит «схлопывание» пленки и счет прекращается.
Рис. 6. Низкочастотная часть спектра изгибных колебаний поверхностных атомов пленки при различных температурах.
Вертикальной пунктирной линией показана минимальная частота утт = 0.0166 ТИг
Появление «заторможенной» моды, относящейся к самой длинной изгибной волне, является обязательным условием для развития процесса разрушения в рассматриваемой модели. В самом деле, в другом направлении вдоль оси у, где тоже выставлены периодические граничные условия, временная зависимость Ьу(() имеет особенность только при временах, когда процесс разрушения в направлении г уже развился, что демонстрирует пунктирная линия на рис. 3. Геометрия рассматриваемой задачи такова, что максимальная длина волны в направлении у в 2.8 раза меньше, чем вдоль оси г. Период ее колебаний, пропорциональный X2 (см. (1) из работы [5]), в восемь раз меньше, чем у аналогичной волны, распространяющейся вдоль направления г и она при всех рассмотренных температурах не становится «заторможенной».
Таким образом, можно выделить три этапа при моделировании методом МД процесса разрушения тонкой пленки:
1. При повышении температуры образуется «заторможенная» изгибная волна, период колебаний которой больше временного интервала существования стационарного состояния (0.6 м).
2. Постоянное со временем увеличение амплитуды этой волны и, соответственно, кривизны пленки до критического значения, когда становится возможным флуктуационное притягивающее взаимодействие между атомами с противоположных склонов гребня волны.
3. Процесс развивается экспоненциально и приводит к разрушению пленки.
ЗАКЛЮЧЕНИЕ
Таким образом, методом молекулярной динамики в работе показана роль низкочастотных изгибных колебаний в процессе разрушения ультратонкой пленки. Установлено, что на начальном этапе разрушения линейное уменьшение размера базового кристаллита обусловлено ростом амплитуды изгибной волны, максимальной для выбранной геометрии задачи. При увеличении амплитуды этой волны растет изгиб пленки до
критической величины, при которой значения элемента тензора внутренних напряжений о^ становятся отрицательными. Это объясняется возникновением флуктуационного притягивающего взаимодействия между атомами с противоположных «рукавов» пленки. После чего процесс развивается лавинообразно и происходит разрушение пленки.
Работа выполнена с использованием кластера «Уран» СКЦ ИММ УрО РАН в рамках темы НИР УдмФИЦ УрО РАН «Теоретические исследования электронных, магнитных, решеточных и транспортных свойств слоистых и наноструктурированных систем» АААА-А17-117022250041-7.
СПИСОК ЛИТЕРАТУРЫ
1. Evans J. W., Thie P. A., Bartelt M. C. Morphological evolution during epitaxial thin film growth: Formation of 2D islands and 3D mounds // Surface Science Reports, 2006, vol. 61, iss. 1-2, pp. 1-128.
2. Venables J. A., Seguin J. L., Suzanne J., Bienfait M. Crystallograghy and growth modes of thick physisorbed films on graphite // Surface Science, 1984, vol. 145, iss. 2-3, pp. 345-363.
3. Venables J. A., Spiller G. D. T., Hanbucken M. Nucleation and growth of thin films // Reports on Progress in Physics, 1984, vol. 47, no. 4, pp. 399-459.
4. Venables J. A. Atomic processes in crystal growth // Surface Science, 1994, vol. 299-300, pp. 798-817.
5. Саламатов Е. И., Долгушева Е. Б. Дисперсия длинноволновых фононов в тонких пленках. Метод молекулярной динамики // Химическая физика и мезоскопия. 2019. T. 21, № 3. С. 420-427.
6. Cross M. C., Lifshits R. Elastic wave transmission at an abrupt junction in a thin plate with application to heat transport and vibrations in mesoscopic systems // Physical Review B, 2001, vol. 64, iss. 8, pp. 085324.
7. Salamatov E. I. Elastic-energy propagation in mesoscopic systems // Physica Status Solidi B, 2013, vol. 250, iss. 9, pp. 1932-1936.
8. Salamatov E. I. Simulation of long-wave phonon (^ > b) scattering at geometric imperfections in nanowires by FDTD method // Physica Status Solidi B, 2012, vol. 249, iss. 1, pp. 113-117.
9. Киселев С. П. Метод молекулярной динамики в механике деформированного твердого тела // Прикладная механика и техническая. 2014. Т. 55, № 3. С. 113-139.
10. Marquardt D. W. An Algorithm for Least-Squares Estimation of Nonlinear Parameters // Journal of the Society for Industrial and Applied Mathematics, 1963, vol. 11, iss. 2, pp. 431-441.
Atomic Mechanism of Thin Film Destruction at High Temperatures. Molecular Dynamics Method
Salamatov E. I., Dolgusheva E. B.
Udmurt Federal Research Center, Ural Branch of the Russian Academy of Science, Izhevsk, Russia
SUMMARY. The atomic mechanism of destruction of an ultrathin aluminum film consisting of five atomic layers is investigated by the method of molecular dynamics, and factors pointing to the beginning of this process are revealed. Special attention is paid to the low-frequency part of the vibrational spectrum of the system. In this connection, a previous symmetry analysis of atomic displacements is performed, which allows independent calculation of the frequency spectra corresponding to the bending and compression vibrations of the film. Thus, it is possible to clearly determine to which branch of the phonon spectrum the considered vibrations with a given frequency V belong, and to find the wavelength X of the corresponding phonon. After determining the temperature range (0 - 600 K) in which the film remains stable for a time of 0.6 ns, the symmetry approach is used to perform a frequency analysis of the vibrational spectrum, starting from a minimum frequency Vmin = 0.0166 THz. Based on the obtained temperature dependences of the low-frequency spectrum, the changes in the film size along the z axis and its shape, a possible scenario for the destruction of thin films in MD simulation is proposed. At the initial stage, with increasing temperature, in the system there arises a "retarded" bending mode corresponding to the wave of maximum length whose vibration period is greater than the time of existence of the stationary state (in our case, more than 0.6 ns). Then a constant in time increase of the amplitude of this mode is observed, which results in a decrease of its
length to a critical value, at which there arises an attractive fluctuation interaction between atoms located on opposite slopes of the wave crest, after which the process develops exponentially, leading to the film destruction. Thus, it is shown that low-frequency bending vibrations corresponding to the wave with the maximum length in the system under consideration play a significant role in the process of film destruction.
KEYWORDS: molecular dynamics method, vibration dispersion, phonon spectrum, Fourier transform, "retarded" mode.
REFERENCES
1. Evans J. W., Thie P. A., Bartelt M. C. Morphological evolution during epitaxial thin film growth: Formation of 2D islands and 3D mounds. Surface Science Reports, 2006, vol. 61, iss. 1-2, pp. 1-128. https://doi.org/10.1016/j.surfrep.2005.08.004
2. Venables J. A., Seguin J. L., Suzanne J., Bienfait M. Crystallograghy and growth modes of thick physisorbed films on graphite. Surface Science, 1984, vol. 145, iss. 2-3, pp. 345-363. https://doi.org/10.1016/0039-6028(84)90087-6
3. Venables J. A., Spiller G. D. T., Hanbucken M. Nucleation and growth of thin films. Reports on Progress in Physics, 1984, vol. 47, no. 4, pp. 399-459. https://doi.org/10.1088/0034-4885/47/4/002
4. Venables J. A. Atomic processes in crystal growth. Surface Science, 1994, vol. 299-300, pp. 798-817. https://doi.org/10.1016/0039-6028(94)90698-X
5. Salamatov E. I., Dolgusheva E. B. Dispersiya dlinnovolnovih fononov v tonkih plenkah. Metod molekulyarnoi dinamiki [Dispersion of long-wave phonons in thin films. molecular dynamics method]. Khimicheskaya fizika i mezoskopiya [Chemical Physics and Mesoscopy], 2019, vol. 21, no. 3, pp. 420-427. https://doi.org/10.15350/17270529.2019.3.44
6. Cross M. C., Lifshits R. Elastic wave transmission at an abrupt junction in a thin plate with application to heat transport and vibrations in mesoscopic systems. Physical Review B, 2001, vol. 64, iss. 8, pp. 085324. https://doi.org/10.1103/PhysRevB.64.085324
7. Salamatov E. I. Elastic-energy propagation in mesoscopic systems. Physica Status Solidi B, 2013, vol. 250, iss. 9, pp. 1932-1936. https://doi.org/10.1002/pssb.201349036
8. Salamatov E. I. Simulation of long-wave phonon (^ > b) scattering at geometric imperfections in nanowires by FDTD method. Physica Status Solidi B, 2012, vol. 249, iss. 1, pp. 113-117. https://doi.org/10.1002/pssb.201147186
9. Kiselev S. P. Method of molecular dynamics in mechanics of deformable solids. Journal of Applied Mechanics and Technical Physics, 2014, vol. 55, no. 3, pp. 470-493. https://doi.org/10.1134/S0021894414030109
10. Marquardt D. W. An Algorithm for Least-Squares Estimation of Nonlinear Parameters. Journal of the Society for Industrial and Applied Mathematics, 1963, vol. 11, iss. 2, pp. 431-441. https://doi.org/10.1137/0111030
Саламатов Евгений Иванович, доктор физико-математических наук, главный научный сотрудник, УдмФИЦ УрО РАН, тел. 8(3412)213966, e-mail: esalama2i@gmail. com
Долгушева Елена Борисовна, кандидат физико-математических наук, старший научный сотрудник, УдмФИЦ УрО РАН, тел. 8(3412)213966, e-mail: elena@udman.ru
