Континуально-атомистическое моделирование абляции и образования наночастиц при воздействии фемтосекундного лазерного импульса на тонкие фольги Текст научной статьи по специальности «Физика»

УДК 544.032.65

В. Б. Фокин, М. Е. Поварницын, П. Р. Левашов КОНТИНУАЛЬНО-АТОМИСТИЧЕСКОЕ МОДЕЛИРОВАНИЕ АБЛЯЦИИ И ОБРАЗОВАНИЯ НАНОЧАСТИЦ ПРИ ВОЗДЕЙСТВИИ ФЕМТОСЕКУНДНОГО ЛАЗЕРНОГО ИМПУЛЬСА

НА ТОНКИЕ ФОЛЬГИ

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

образование наночастиц.

В работе рассмотрена гибридная континуально-атомистическая модель для моделирования лазерной абляции в металлах. Движение ионов описывается методом молекулярной динамики, для электронов решается одномерное уравнение энергии, при этом учитывается обмен энергией между этими двумя подсистемами. Численно определен порог абляции и обнаружено возрастание глубины абляционного кратера при повышении интегральной плотности потока энергии лазерного излучения в разумном соответствии с экспериментальными данными. Также продемонстрировано образование наночастиц.

Keywords: femtosecondlaserpulses, ablation, moleculardynamics, two-temperature model, formation of nanoparticles.

In this work we consider the hybrid continuum-atomistic model for the simulation of laser ablation in metals. The motion of ions is governed by the molecular dynamics method, for the electrons we solve the one-dimensional energy equation, wherein the energy exchange between these two subsystem is taken into account. The ablation threshold is determined numerically as well as the rise of the ablation crater depth at increasing laser fluence is revealed in reasonable agreement with experimental data. Also wedemonstrate the formation of nanoparticles.

Лазерная обработка материалов широко применяется в современных технологических процессах [1, Ошибка! Неизвестный аргумент ключа., Ошибка! Неизвестный аргумент ключа.],

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

Для моделирования взаимодействия лазерного излучения с металлами в основном используются два подхода: двухтемпературная гидродинамическая модель (ГД)[4,5,6] и метод молекулярной динамики (МД) [7]. У каждого из этих методов есть свои преимущества и недостатки. Так, в ГД необходимо использование дополнительных моделей, описывающих кинетику фазовых превращений, фрагментации вещества, эффекты, связанные с поверхностным натяжением; следует отметить, что многие из этих свойств недостаточно изучены. Кроме того, гидродинамический подход неприменим на нанометровых масштабах. С другой стороны, метод МД для своей работы требует лишь знания потенциала межчастичного взаимодействия, а потому позволяет естественным образом учитывать структуру кристаллической решетки, влияние примесей и наличие дислокаций, кинетику фазовых переходов, явления нуклеации и откола при растяжении вещества и многие другие процессы. При этом адекватность их учета в веществе полностью определяется качеством этого потенциала. Важно отметить, что область применимости многих из существующих потенциалов весьма ограничена. Кроме того, в МД

чрезвычайно сложно учесть квантовую природу электронов в металлах.

В данной работе используется гибридная континуально-атомистическая модель, которая объединяет достоинства ГД и МД методов и основана на подходе работы [8]. Такая комбинированная модель воспроизводит динамику вещества, как во время лазерного воздействия, так и на более поздних временах, учитывая поглощение излучения, различные температуры электронов и ионов, эффекты, связанные с теплопроводностью, процессы плавления решетки, нуклеацию и откол вещества.

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

Быстрый лазерный нагрев электронов в зоне проводимости приводит к появлению неравновесного состояния с «горячими» электронами и «холодной» решеткой [8]. Более того, неравновесное тепловое состояние характеризуется резким градиентом электронной температуры, что является причиной появления силы электронного давления [10], существенной, по крайней мере, на ранних стадиях лазерного нагрева (порядка нескольких пикосекунд). Таким образом, для каждого атома в системе мы решаем уравнение движенияв следующем виде:

d 2r

_ F, vT + Ff .

(1)

d(peei) , d(peeiu) ^ _ du

' С2

Здесь т,— масса / -го атома, г,— радиус-вектор атома, Р,— сила, возникающая вследствие взаимодействия /-го атома с окружающими атомами, —коэффициент обмена, возникающий вследствие отсутствия теплового равновесия между электронами и ионами [8], — сила электронного давления, \,Т— тепловая скорость/ -го атома (скорость атома относительно центра масс элементарного объема V, содержащего этот атом). Объем V достаточно мал по сравнению с размером мишени, но,тем не менее, содержит много атомов N >> 1. Тепловая скорость атома/ может быть выражена следующим образом:

N

T ^ c

V _f v i - v .

i_1

(2)

атома,

Здесь ={ы1, V 1, w¡} — скорость Vе = |мс,Vе, wе |— скорость центра масс системы

атомов, заключенных внутри объема V, которую можно вычислить следующим образом:

Vе =£(т,.V,.)/¿т,. . (3)

,=1 ,=1

Ранее в работе [8] было получено выражение для коэффициента который описывает

температурную релаксацию электронной и ионной подсистем. Вначале вводится температура атомов внутри объема Vкак средняя кинетическая энергия, приходящаяся на атом, следующим образом:

То"_ Ü^N f m (vT )2

(4)

где kB—постоянная Больцмана. Затемкоэффициент обмена вычисляется по формуле

_yV ( - Ton) (5)

bi — bV N _ „ , , rrr, (5)

f m (vT )2

3kBNTion

и полагается одинаковым для всех атомов, принадлежащих объему V [8].

Наконец, сила давления электронного газаследует из градиента электронного давления, и

вычисляется как = т УР1. В данной работе мы

Р

пренебрегаем поперечными эффектами

теплопроводности и электронного давления, и градиент давления имеет только одну г-компоненту, нормальную к свободной поверхности фольги. Надо отметить, что модели [8, 11] не гарантируют сохранения полной энергии системы, так как уравнение теплопроводности для свободных электронов не учитывает движение вещества и работу электронного давления. Для решения этих проблем мы описываем эволюцию подсистемы свободных электронов с помощью уравнения энергии. Это уравнение записывается в форме Эйлера следующим образом:

dt

dz

dz

- To-)+ql (t, z)

МД потенциал

Для МД-моделирования алюминия, используется EAM-потенциал Жаховского [12]. Потенциал разработан на основе

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

В качестве уравнения состояния электронов используетсятепловая часть конечнотемпературной модели Томаса-Ферми [13], которая дает соотношения в форме Pei(p, Te), eel(p, Te), используемые для замыкания уравнения (6).

Результаты моделирования и их обсуждение

В нашей работе моделировался титан-сапфировый лазер с длиной волны 800 нм и длительностью импульса 100 фс, и результаты моделирования сравнивались с

экспериментальными данными [14]. Временной профиль интенсивности импульса имел Гауссову форму I(t) = ^expt-ln^Xt-ta)2/?2^ пиковой интенсивностью I0. Соотношение между интегральной плотностью падающего лазерного излучения и пиковой интенсивностью выглядит

следующим образом: Finc _ I0r^/ln(16) .

Моделирование начинается в момент времени 0 пс, максимум интенсивности приходится на момент времени t0 = 2 пс.

В начальный момент времени свободная поверхность мишени расположена в точке с координатой z = 0, а импульс воздействует на мишень слева. Толщины мишени должно быть достаточно, чтобы избежать воздействия волн давления, отраженных от правой части мишени, на развитие процесса абляции. Для МД-моделирования это соответствует толщине 2 мкм, а поперечное сечение фольги составляло примерно 8 х 8 нм (20 периодов кристаллической решетки типа ГЦК). Указанной толщины мишени достаточно, чтобы исследовать начальную стадию лазерной абляции вплоть до времен 500 пс.

Было проведено моделирование взаимодействия лазер - металл для случая интегральной плотности падающего излучения, близкой к порогу абляции (для алюминия порог абляции около 0.2 Дж/см2 [14]). Для случая интегральной плотности падающего излучения Fmc = 0.5 Дж/см2 (Fabs = 0.065 Дж/см2) на рис. 1 представлена x-t диаграмма плотности вещества.

Полосатая структура плотности на диаграмме в моделировании возникает из-за того,

i-1

что кристаллическая решетка ГЦК ориентирована вдоль оси г. Отчетливо видно, что формируется волна сжатия, которая движется со скоростью около 5.5 км/с, это значение примерно соответствует скорости звука в алюминии при нормальных условиях. Расширение вещества из области высокого давления в вакуум приводит к формированию волны растяжения и фрагментации мишени на несколько слоев в момент времени 20 пс. Этот слой отлетает от мишени со скоростью около 0.5 км/с. Откол слоев происходит в жидкой фазе и, таким образом, вклад в формирование кратера вносит нуклеация в разреженной метастабильной жидкой фазе при отрицательном давлении.

Плотность, г/см3

3.5

140

120

100

о

[=

□Г 80

Ф О. 60

ш

40

20

0

■ 1 Ън В 1

1 11 | I ! У7

I Si ■ i ¡11 \ ■

PI

ШИВ

\ 1 У .s * / ^^^

I

3

2.5

2

1.5 1

0.5 0

500

-200 -100 0 100 200 300 400 Координата по оси Ъ, нм Рис. 1 - х^ диаграмма плотности при воздействии на фольгу лазерного импульса с интегральной плотностью излучения Ж = 0.5 Дж/см2. Пунктиром показан фронт плавления, сплошная стрелка - направление распространения ударной волны, штрих-пунктирная стрелка -направление распространения волны разрежения

В другом расчетеинтегральная плотность падающего излучения была увеличена до = 2.0 Дж/см2 = 0.4 Дж/см2). Передние слои

отколовшегося материала находятся в состоянии жидкость-газ (рис. 2).

200

фазовом

150

100

50

-600 -400 -200 0 200 Координата по оси Z, нм

400

Рис. 2 - х-^диаграмма плотности при воздействии на фольгу лазерного импульса с интегральной плотностью излученияЖ = 2.0Дж/см2

На этой диаграмме явно различимы отдельные наночастицы. Более глубокие слои разрушаются в жидкой фазе. Общая глубина абляции для данной интенсивности составляет около 100 нм.

Далее, мы промоделировали процесс фрагментации тонкой пенки толщиной 100 нм под воздействием излучения с интегральной плотностью излучения Жтс = 1 Дж/см2. Фрагмент системы атомов в момент времени 450 пс выглядит следующим образом (рис.3).

Рис. 3 - Фрагмент системы атомов в момент времени 450 пс, интегральная плотность излучения Ж = 1.0 Дж/см2

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

Наконец, мы нанесли полученные в данной работе точки на график зависимости глубины кратера абляции от интенсивности, приведенный в работе [14]. Глубину кратера О, получающуюся в ходе такого эксперимента, можно найти, подсчитав при моделировании суммарную массу, которая отделилась от образца в процессе фрагментации:

т

О =-, (7)

р0 £

где т— значение аблированной массы, полученной в расчетах, р0— начальная плотность алюминия, £— площадь поперечного сечения образца, используемого в моделировании. Данная зависимость приведена на рис. 4, точки, полученные в данной работе, выделены красным цветом.

Интегральная плотность излучения лазера F, Дж/см

Рис. 4 - Глубина абляции алюминия D в зависимости от интегральной плотности падающего излучения Finc для фемтосекундных импульсов на длине волны 800 нм. Квадратики

— моделирование данной работы, треугольники

— экспериментальные данные [14]

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

Выводы

В данной работе представленагибридная двухтемпературная модель для расчета лазерной абляции металлов. Поглощение лазерной энергии, электронная теплопроводность и электрон-ионное взаимодействие учитывается с применением широкодиапазонных моделей оптических и транспортных свойств. Для моделирования сил межатомного взаимодействия используется потенциал в форме ЕАМ. Показано, что наш подход дает разумное соответствие с экспериментальными данными для интегральных плотностей энергии лазера 0.5 и 2.0 Дж/см2.Разработанный метод требует значительных вычислительных ресурсов, но позволяет воспроизводить процессы в веществе на атомарном уровне. Среди них такие, как нуклеация в метастабильной жидкой фазе,формирование газожидкостной смеси, образование наночастиц, для которых поверхностные эффекты играют определяющую роль. Данный подход позволяет улучшить пониманиеявлений лазерной абляции и может быть использован для построения кинетических моделей разрушения в гидродинамическом коде в будущем.

Работа выполнена при финансовой

поддержке РФФИ, гранты № 13-02-91057, 13-0801179.

Литература

1. А. Г. Григорьянц, И. Н. Шиганов, А. И. Мисюров. Технологические процессы лазерной обработки: учебное пособие для ВУЗов. Изд-во МГТУ им. Н.Э. Баумана, Москва, 2006. 664 с.

2. М. Ю. Зайцева, Л. Н. Абуталипова.Вестник Казан. технол. ун-та, 6,93-94 (2013)

3. А. А. Ефремова, Р. М. Гарипов.Вестник Казан. технол. ун-та, 3, 85-87 (2012)

4. J. P. Colombier, P. Combis, F. Bonneau, R. Le Harzic, and E. Audouard, Phys. Rev. B.,71, 165406 (2005)

5. S. Anisimov, N. Inogamov, Y. Petrov, V. Khokhlov, V. Zhakhovskii, K. Nishihara, M. Agranat, S. Ashitkov, and P. Komarov, Applied Physics A,92, 939 (2008)

6. M. E. Povarnitsyn, N. E. Andreev, E. M. Apfelbaum, T. E. Itina, K. V. Khishchenko, O. F. Kostenko, P. R. Levashov, and M. E. Veysman, Applied Surface Science,258, 9480 (2012)

7. Y. Cherednikov, N. A. Inogamov, and H. M. Urbassek, Phys. Rev. B,88, 134109 (2013)

8. D. S. Ivanov and L. V. Zhigilei, Phys. Rev. B,68, 064114 (2003)

9. M. S. Daw and M. I. Baskes, Phys. Rev. B,29, 6443 (1984)

10. J. Chen, D. Tzou, and J. Beraun, International Journal of Heat and Mass Transfer,49, 307 (2006)

11. Y. Gan and J. K. Chen, Applied Physics Letters,94, 201116(2009)

12. V. Zhakhovskii, N. Inogamov, Y. Petrov, S. Ashitkov, andK. Nishihara, Applied Surface Science,255, 9592 (2009)

13. O. P. Shemyakin, P. R. Levashov, L. R. Obruchkova, and K. V. Khishchenko, J. Phys. A: Math. Theor.,43, 335003 (2010)

14. M. Hashida, A. Semerok, O. Goberta, G. Petiteb, and J.-F. Wagner, in SPIE Proceedings, 4423 pp. 178-185 (2001)

© В. Б. Фокин - м.н.с., ОИВТ РАН, Москва, vladimirfokin@mail.ru; М. Е. Поварницын - с.н.с., ОИВТ РАН, Москва, povar@ihed.ras.ru; П. Р. Левашов - с.н.с., ОИВТ РАН, Москва, pasha@ihed.ras.ru.

© V. B. Fokin - junior researcher, JIHT RAS, Moscow, vladimirfokin@mail.ru; M. E. Povarnitsyn - senior researcher, JIHT RAS, Moscow, povar@ihed.ras.ru; P. R. Levashov - head of laboratory, JIHT RAS, Moscow, pasha@ihed.ras.ru.