МОДЕЛИРОВАНИЕ ВОЗДЕЙСТВИЯ НА СТЕКЛО ПОСЛЕДОВАТЕЛЬНОСТИ ЛАЗЕРНЫХ ИМПУЛЬСОВ С РАЗНЫМИ ДЛИНАМИ ВОЛН Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 22, № 6, 2017

Моделирование воздействия на стекло последовательности лазерных импульсов с разными длинами волн

В. П. Жуков1'2'*, Н.М. Булгакова3'4, М. П. Федорук1'5 1Институт вычислительных технологий СО РАН, Новосибирск, Россия 2Новосибирский государственный технический университет, Россия 3Центр Хайлэйс Института физики АН ЧР, Прага, Чешская республика 4Институт теплофизики СО РАН, Новосибирск, Россия 5 Новосибирский государственный университет, Россия *Контактный e-mail: zukov@ict.nsc.ru

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

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

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

Введение

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

© ИВТ СО РАН, 2017

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

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

Ввиду сложности и высокой стоимости натурного эксперимента возникает необходимость численного моделирования предложенной схемы обработки материала дихроматическими импульсами. Описанные в настоящей работе результаты моделирования позволяют количественно оценить ожидаемый эффект и предложить примерные параметры лазерных импульсов. Исследуется схема, в которой второй импульс имеет такую задержку, что плазма электронов проводимости, созданная первым импульсом, еще не успевает рекомбинировать, но импульсы при этом не перекрываются. В качестве материала выбран плавленый кварц с шириной запрещенной зоны Ед0 = 9 эВ и плотностью атомов р0 = 6.6 • 1022 см-3.

1. Математическая модель и метод решения

Для моделирования описанной выше схемы воздействия последовательности импульсов на стекло использовались следующие уравнения [2-5]:

1Ж „8^ Е

-с-Ш — Б = -^ +гоШ - Vе!Е|2, (1)

Б = п2Е + -^П2П2\Е\2Е, (2)

4тт

е V

гш\ = — Е +--, J = — реV, е > 0, (3)

те тс

^ = WpJ + — р/Пг, (4)

w =_е2 Тср (р о — р)\Щ2_

2теЕдо(1 + т/те)(1 + —2г2)Ро (1 + |Е|2/(4Е2)), ( )

Ед0т—2 1д Н — ^

Е2 = д02 , 1— — г-Н = — тЬЕ. е2 с оъ с

Здесь все физические поля представлены с помощью комплексных амплитуд. Например, реальное физическое электрическое поле описывается как = (Ее-гш4+Е*ешЬ)/2, где — — частота несущей волны. Среднее по периоду значение квадрата поля равно (Е2^) = \Е\2/2, где \Е\2 = (ЕЕ*) — квадрат модуля комплексной амплитуды.

Остальные обозначения следующие: V и р — гидродинамическая скорость и плотность электронов в зоне проводимости; те — масса электрона в зоне проводимости; которая полагается равной массе электрона в вакууме; т = 0.5те — приведенная масса; п — показатель преломления; п2 — коэффициент в эффекте Керра; тс = 1.28 фс — время столкновения электрона с атомами, = 150 фс — время захвата электронов зоны проводимости в ловушки; — скорость ударной ионизации.

Для скорости фотоионизации использовалось выражение Шр1 = Шр10Шк(ш, (1 — р/ро), где Шк — зависящая от ш и 1Е|2 скорость фотоионизации, рассчитанная по теории Келдыша [6]. Величина Шк вычислялась численно без использования метода перевала, как в [6]. Параметр Шр20 в настоящей работе подбирался таким, чтобы обеспечить соответствие предыдущим работам [2-5]. В наших расчетах Шр^ = 3.6^ 1033 с-1 •см-3 при длине волны излучения 800 нм, р = 0 и 1Е |2 = Е2.

Энергия, поглощаемая из излучения при фотоионизации С, рассчитывалась следующим образом. Согласно [6] Шк = ып, где и>га — скорость ионизации в и-фотонном

п

процессе. Соответственно, С = ЭДр/0(1 — р/р0) ^ пНштп.

п

Приведенные уравнения решались в двумерной аксиально-симметричной (г — г)-гео-метрии для линейно-поляризованного импульса с помощью конечно-разностной схемы на сдвинутых сетках, подобной схеме, описанной в [3-5] для несколько других уравнений. В частности, аппроксимация пространственных производных аналогична схеме Ии [7] и совпадала с [3-5]. Аппроксимация уравнений (1)-(3) по времени имеет вид

1 вга+1 — Бга ш Бга+1 + Бга vn + vra+1 „ + тт„ 8^ Еп ---г----=----ерп+1 + гсШга--С

. , с 2 с 2 с |Era|2

п

Dn = n2En + —n2^2|Era-1|2Era

е vn

iuvn = — Era + —, me тс

Подчеркнем, что при возникающих в данной задаче величинах плотности р такая неявная схема необходима. Расчет по явным схемам практически невозможен из-за слишком малого шага по времени т, который необходим в связи как с устойчивостью, так и с точностью расчета. Как видно, схема имеет первый порядок аппроксимации по т. Уравнение для плотности также решалось простейшим методом Эйлера первого порядка по времени. Пространственные производные аппроксимировались со вторым порядком.

Оценки устойчивости схемы в предположении постоянства р в случае декартовых координат показывают, что схема устойчива при условии Куранта ст < const h (h — наименьший пространственный шаг сетки, const — коэффициент порядка 1). Расчеты подтверждают это условие. Расчеты также показывают, что при выполнении условия устойчивости шаг по времени не оказывает существенного влияния на точность расчета, которая определяется в основном шагом пространственной сетки. Заметим, что, естественно, наибольшая погрешность возникает при расчете импульса с меньшей длиной волны. При этом использовалась одна и та же сетка как для коротковолновых, так и для длинноволновых импульсов. В расчетах шаг пространственной сетки по оси z составлял 0.025 мкм, а по радиусу — 0.008 мкм. Длина волны коротковолнового излучения лазера равнялась 0.4 мкм. Погрешность расчетов на такой сетке оценивается в 5%.

Начальные условия и вид ж-компоненты электрического поля (линейная поляризация), задаваемой на границе z = 0, описаны ниже. На внешних границах использовались неотражающие (поглощающие) граничные условия. Существует множество вариантов таких условий [8-10]. В нашем случае интерес представляют только процессы, происходящие в фокусе. Кроме того, нас интересуют короткие импульсы. Отраженное от границ излучение возвращается в фокусную область после того, как все интересующие нас процессы произошли. Соответственно, оно не оказывает большого влияния на исследуемые явления. Поэтому в настоящей работе реализован простейший вариант неотражающих граничных условий. А именно, после каждого шага по времени электрическое и магнитное поля, токи, скорости свободных электронов в точках, удаленных от границы на расстояние А, меньшее некоторого числа dp^M ~ 1 мкм, умножались на число (1 + тд0(1 — А/dpLM))-1. Здесь т — шаг по времени, фс. Типичные значения параметров таковы: dpLM = 1.5 мкм, д0 = 5. Они, а также параметры формы расчетной области, радиус которой уменьшался по мере приближения к фокусу, выбирались экспериментально из расчетов. Отметим, что радиус узкой части расчетной области составлял 20-25 мкм при радиусе пучка в фокусе (перетяжке) в 1 мкм. В фокусе концентрируются волны, приходящие от областей пространства, удаленных по радиусу. Поэтому значительное уменьшение наименьшего радиуса расчетной области не представляется возможным. Радиус расчетной области при z = 0 определяется шириной задаваемого пучка w и должен быть не меньше 2-3 w.

2. Результаты расчетов

Рассмотрим результаты следующего расчета. В начальный момент времени t = — t0, t0 = 500 фс электромагнитное поле и электроны в зоне проводимости отсутствуют. Первый лазерный импульс с длиной волны в вакууме Л = 400 нм (п = 1.47, п2 = 0.35 в случае плавленого кварца [11]) инициируется условием для ж-компоненты электрического поля на границе = 0 вида

Ех = Е0 e-t2/tl >.

Здесь k = 2пп/Л — волновой вектор излучения в стекле; f ~ 120 мкм приблизительно соответствует фокусному расстоянию; w = 20.4 мкм — радиус лазерного пучка на границе. Перетяжка для пучка с такими параметрами в линейном режиме (Е0 ^ 0) равна Wf = 0.509 мкм. Длительность импульса tb = 140 фс. Амплитуда импульса Е0 соответствовала энергии 0.006 мкДж.

Имеется обычная картина распространения импульса: излучение лазера сходится к фокусу, и в окрестности последнего (в результате фотоионизации) появляются электроны зоны проводимости. Они поглощают падающее излучение и сами производят новые электроны благодаря ударной ионизации. Параметры расчета соответствуют слабому импульсу, при котором максимальная поглощенная энергия Еаьт = max Еаь ~

r,z

608 Дж/см3, где

Е0ь = / (j*^ + G) d,

-to

Этой энергии недостаточно для какой-либо остаточной модификации стекла. Распределение Еаь показано на рис. 1. Аналогичное распределение имеет плотность плазмы р.

При этом тахр ^ 0.0029р0. Отметим, что в этом случае вклад фотоионизации в плотность плазмы более чем в два раза превосходит вклад ионизации электронным ударом.

При Ь > 800 фс импульс покидает фокусную область, напряженность электрического поля становится малой. Взаимодействие излучения с материалом прекращается. Плотность плазмы постепенно падает за счет рекомбинации (последний член в (4)).

В момент времени Ь = 800 фс расчет импульса с Л = 400 нм прекращается и начинается расчет импульса с Л = 800 нм. Для этой длины волны п = 1.4533, п2 = 0.248.

О 0.2 0.4 Г, мкм О 0.2 0.4 0.6 г, мкм

О 1 Г, мкм 0 1 2 3 4 Г, мкм

Рис. 1. Распределение поглощенной энергии (Дж/см3) после прохождения первого импульса с Л = 400нм (а); двух последовательных импульсов с Л = 400 и 800нм (б); одиночного импульса с Л = 800 нм и энергией 0.12 мкДж (в) и одиночного импульса с Л = 800 нм и энергией 1.5 мкДж (г)

Второй импульс имеет параметры = 180 фс, = 10.7мкм, / ^ 120 мкм. Соответственно перетяжка Wf = 2 мкм. Таким образом, второй импульс намного шире первого. Поэтому проблем с совмещением фокусных областей первого и второго импульсов в реальном эксперименте не должно возникать. Энергия этого импульса равна 0.12 мкДж. Реально расчет второго импульса начинается раньше. В момент времени Ь = 800 фс интенсивность второго импульса на границе г = 0 достигает максимума. Но в начальные моменты времени второй импульс находится далеко от фокусной области и не взаимодействует с материалом.

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

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

0.002

0.008

0.006

0.004

О

300

800

1300

и фс

Рис. 2. Зависимость максимальной в пространстве плотности плазмы от времени в результате распространения последовательности из двух импульсов

Рис. 3. Эволюция распределения (интенсивности излучения) при распространении

второго импульса

щенной в результате прохождения двух импульсов энергии существенно выше, чем каждого импульса в отдельности (рис. 1), и ее достаточно для возникновения остаточной модификации стекла. При этом размер области модификации мал и определяется коротковолновым импульсом. Заметим, что в случае одиночного импульса с Л = 800 нм и теми же параметрами плотность поглощенной энергии очень мала, а размер области поглощения значительно больше, чем в случае последовательности двух импульсов. Увеличение энергии одиночного импульса с Л = 800 нм приводит к значительному увеличению размера области поглощения. Так, характерный размер области поглощенной энергии (см. рис. 1) в случае импульса с энергией 1.5 мкДж в разы превосходит этот размер в случае двойного импульса. При этом поглощенная энергия по-прежнему далека от порога остаточной модификации.

Заключение

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

Благодарности. Исследование выполнено при финансовой поддержке РФФИ (грант № 15-01-02432).

Список литературы / References

[1] Potemkin, F.V., Bravy, B.G., Bezsudnova, Yu.I., Mareev, E.I., Starostin, V.M., Platonenko, V.T., Gordienko, V.M. Overcritical plasma ignition and diagnostics from oncoming interaction of two color low energy tightly focused femtosecond laser pulses inside fusied silica // Laser Phys. Lett. 2016. Vol. 13, No. 4. Paper 045402c. 8 p.

[2] Couairon, A., Sudrie, L., Franco, M. et al. Filamentation and damage in fused silica induced by tightly focused femtosecond laser pulses // Phys. Review B. 2005. Vol. 71, No. 12. Paper 125435. 11 p.

[3] Bulgakova, N.M., Zhukov, V.P., Meshcheryakov, Yu.P., Gemini, L., Brajer, J., Rostohar, D., Mocek, T. Pulsed laser modification of transparent dielectrics: What can be foreseen and predicted in numerical experiments // J. of the Optical Soc. of America. B. 2014. Vol. 31, No. 11. P. C8-C14.

[4] Zhukov, V.P., Rubenchik, A.M., Fedoruk, M.P., Bulgakova, N.M. Interaction of doughnut-shaped laser pulses with glasses // JOSA B. 2017. Vol. 34, No. 2. P. 463-471.

[5] Булгакова Н.М., Жуков В.П., Федорук М.П. Численное моделирование распространения фемтосекундного лазерного импульса в нелинейных средах // Вычисл. технологии. 2012. Т. 17, № 4. С. 14-28.

Bulgakova, N.M., Zhukov, V.P., Fedoruk, M.P. Numerical modeling of propagation of femtosecond laser pulses in non-linear media // Comput. Technologies. 2012. Vol. 17, No. 4. P. 14-28. (In Russ.)

[6] Келдыш Л.В. Ионизация в поле сильной электромагнитной волны // ЖЭТФ. 1964. Т. 47, вып. 5(11). С. 1945-1957.

Keldysh, L.V. Ionization in the field of a strong electromagnetic wave // Soviet Phys. JETP. 1965. Vol. 20, No. 5. P. 1307-1314.

[7] Yee, K.S. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media // IEEE Trans. on Antennas and Propagation. 1966. Vol. 17. P. 585-589.

[8] Berenger, J.-P. A perfectly matched layer for the absorption of electromagnetic waves // J. of Comput. Phys. 1994. Vol. 114. P. 185-200.

[9] Petropoulos, P.G. Reflectionless sponge layers as absorbing boundary conditions for the numerical solution of Maxwell equations in rectangular, cylindrical, and spherical coordinates // SIAM J. Appl. Math. 2000. Vol. 60, No. 3. P. 1037-1058.

[10] Deinega, A., Valuev, I. Long-time behavior of PML absorbing boundaries for layered periodic structures // Comput. Phys. Commun. 2011. Vol. 182. P. 149-151.

[11] Milam, D. Review and assessment of measured values of the nonlinear refractive-index coefficient of fused silica // Appl. Optics. 1998. Vol. 37, No. 3. P. 546-550.

[12] Grojo, D., Gertsvolf, M., Lei, S., Barillot, T., Rayner, D.M., Corkum, P.B.

Exciton-seeded multiphoton ionization in bulk SiO2 // Phys. Review B. 2010. Vol. 81. Paper 212301. 4 p.

[13] Peng, J., Grojo, D., Rayner, D.M., Corkum, P.B. Control of energy in femtosecond laser dielectric interactions // Appl. Phys. Lett. 2013. Vol. 102. Paper 16105. 3 p.

Поступила в 'редакцию 9 ноября 2017 г.

56

B. n. }KyKOB, H. M. BynraKOBa, M. n. QegopyK

Modelling of impact of a sequence of laser pulses with different wavelengths on glass

Zhukov, Vladimir P.1,2'*, Bulgakova, Nadezhda M.3,4, Fedoruk, Mikhail P.1'5

institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia 2Novosibirsk State Technical University, Novosibirsk, 630073, Russia 3HiLASE Centre, Institute of Physics ASCR, Prague, 25241, Czech Republic 4Kutateladze Institute of Thermophysics, SB RAS, Novosibirsk, 630090, Russia 5Novosibirsk State University, Novosibirsk, 630090, Russia * Corresponding author: Zhukov, Vladimir P., e-mail: zukov@ict.nsc.ru

The following scheme of materials fabrication by a sequence of two laser pulses with different wavelength is investigated by numerical modelling. At first, the glass is irradiated by a shortwave pulse, which produces seed electrons in the conductive zone of a small size (in order to several cubic wavelengths) volume due to photoionization. The long wavelength second pulse provides the further increase of conductive electrons density due to impact ionization and high absorption of radiative energy. The absorbed energy is so large, so it provides the modification of material.

Thereby the material modification is created in a small volume of space, which is determined by the smallest wavelength. The energies of both pulses could be well below the energy required for the modification by single pulse. The modelling was done in the frame of nonlinear Maxwell equations supplemented by hydrodynamic type equations for the plasma of conductive zone electrons. The features of used finite-difference scheme are also discussed.

Keywords: femtosecond laser pulses, dichromatic laser pulse, nonlinear Maxwell's equations.

Acknowledgements. The work is financially supported by the RFBR (grant No. 1501-02432).

Received 9 November 2017

© ICT SB RAS, 2017