Численное моделирование распространения фемтосекундных лазерных импульсов в нелинейных средах Текст научной статьи по специальности «Физика»

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

Том 17, № 4, 2012

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

u >к

в нелинейных средах*

Н.М. Булгакова1, В. П. Жуков2, М.П. Федорук2 1 Институт теплофизики СО РАН, Новосибирск, Россия 2Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: nbul@itp.nsc.ru, zukov@ict.nsc.ru, mife@ict.nsc.ru

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

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

Введение

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

Для моделирования распространения ультракоротких лазерных импульсов в нелинейных диспергирующих средах (газах и в ряде случаев в прозрачных твёрдых матери-

* Работа выполнена при финансовой поддержке РФФИ (грант № 12-01-00510-а), Министерства образования и науки РФ и Междисциплинарного интеграционного проекта СО РАН № 68.

алах) обычно используется приближение нелинейного уравнения Шрёдингера (НУШ). Это приближение хорошо зарекомендовало себя в применении к моделированию распространения лазерного излучения в волоконной оптике и в газах, где либо плазма отсутствует (в линиях связи это вредный эффект), либо ее концентрация мала. Для условий модификации оптических материалов фемтосекундными импульсами лазерного излучения длительностью 100 фемтосекунд и менее и энергией выше 0.1 мкДж используется НУШ, дополненное уравнением для плотности свободных электронов [1-7]. Однако исследования в данном направлении далеко не полны. Кроме того, применимость НУШ сильно ограничена, поскольку при выводе этого уравнения используются следующие предположения [8]:

1) малость длины несущей волны по сравнению с характерным пространственным масштабом изменения функции огибающей волнового пакета;

2) отсутствие рассеяния волн на большие углы;

3) малость дивергенции электрического поля.

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

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

Для решения уравнений Максвелла в различных физических условиях разработано множество различных методов (конечных разностей, объёмов, элементов, разложение в ряды Фурье и пр.). Заметим, что в работах других авторов, использующих уравнения Максвелла (см., например, [9-14]), моделируются главным образом явления в волноводах. В этих работах не принимаются во внимание ряд явлений, важных для рассматриваемого здесь случая. Так, не учитывается или учитывается не в полной мере дисперсия показателя преломления, отсутствуют нелинейные эффекты, не рассматриваются процессы ионизации материала. Недавно появилась первая значительная работа с попыткой комбинирования уравнений Максвелла и динамики плазмы в приложении к задачам модификации материалов [15]. Однако предложенная в ней модель не учитывает процессов столкновительной ионизации, которые играют важнейшую роль в формировании плазмы при ультракоротких лазерных импульсах вследствие сильного нагрева газа свободных электронов (до нескольких десятков электронвольт). Учет данного эффекта значительно меняет морфологию лазерно-индуцированной плазмы. В настоящей работе предпринята попытка учёта всего разнообразия процессов, возника-

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

1. Исходные уравнения

Будем описывать распространение лазерного импульса уравнениями Максвелла с учётом частотной дисперсии линейной части вектора поляризации и нелинейной составляющей (эффект Керра). Дополним их уравнениями гидродинамического типа для электронов. Следуя [1, 2, 4, 8], представим электрическое поле в комплексном виде E = ECе-ш1 + ECеш1. Здесь ш — несущая частота волны, EC = EC(r,t) — медленная огибающая. Звёздочкой обозначено комплексное сопряжение. Аналогичные выражения примем для магнитного поля H, тока j, скорости электронов v, электрической индукции D и умноженной на 4п поляризации среды, связанной с m-м осциллятором Pm. Тогда усреднённый по периоду квадрат электрического поля (E2) = |EC |/2. Модуль вектора Пойнтинга, усреднённый по периоду, имеет вид

I = — |(E х H) = А ec х HC + EC х Hc.

4п|( ) 8n 2

В среде с показателем преломления n в случае плоской волны в линейном приближении (|HC| = n|EC|) это даёт I = cn\E2C|/8п. Здесь n — показатель преломления при частоте ш.

Ниже индекс C для краткости будем опускать. Тогда имеем 1 d D .ы 4п. тт 8пТ1Г ^ E

cä-- гш D = - inj +rot H - TW»E' |Ei ■ (1)

D = E + ^ Pm + P„,, (2)

m

- ^Pm = Vm, (3)

m

д V„

dt

dt

- - lшVm = -Ш2т(Pm - BmE) , (4)

c

Pn = 4^n2n2 ((1 - fr) |E|2 + frRaq) E, (5)

Ro=(T?+^)(6)

q = + А (7)

Td Ts

P = IE|2 - p - (8)

Td Ts

d (pv) e ^ pv u J гшpv =--pE - —, (9)

dt me tc

dP = WpJ + WCT - (10)

dt Ttr

j = -pev, e > 0, (11)

1 dH - ¿ШH = -rot E. (12)

c dt - c -

Здесь р — плотность свободных электронов (см-3); Шр/ и Wa — скорости многофотонной и столкновительной ионизации; тс и т4г — характерные времена рассеяния электронов и их рекомбинации; п и п2 — линейный и нелинейный показатели преломления; /г — коэффициент задержки оптического отклика; т3 и та — характерные времена задержки оптического отклика; Ед — эффективный потенциал ионизации; Ед0 — ширина запрещённой зоны. Последнее слагаемое в (1) описывает потери энергии излучения на многофотонную ионизацию.

Для монохроматической волны в линейном приближении уравнения (1)-(5) дают показатель преломления, выражение для которого имеет вид

п2 _ 1 + ^ = 1 + Бт^2

п _1 + т ^ - -1 + т л2 - лт'

т т т т

где Л _ 2пс/ш — длина волны в вакууме.

Уравнения (5)-(8) при _ 0) _ _ 0) _ 0 дают для нелинейной поляризации формулу

Р ш _4п'п2'П2 ( (1 - / )|Е|2 + /г По I |Е2|е-т/т^ 8т(т/т^т ) Е.

Уравнения (3)-(8) позволяют обойти необходимость использования спектральных методов решения исследуемой задачи.

Формулы для Шр1, Ша Ед имеют следующий вид:

Шр1 _ ШРЮ(|Е|2/Е,2)а , (13)

ро

Ш __е2тср(ро - р)|Е21__(14)

2теЕдо(1 + ш/ше)(1+ ^2т2)ро (1 + |Е2|/(4Е?))' ( )

Ед _ Едо (1 + |Е2|/(4Е2)) , (15)

Е? _ ^, (16)

е2

где т — эффективная масса электрона в твердом теле.

Подчеркнем, что значения всех коэффициентов таковы [1, 2, 4], что под |Е2| понимается именно квадрат модуля комплексной амплитуды, а не усреднённый по периоду волны квадрат электрического поля.

Из (1)-(8), (11), (12) можно получить уравнение для плотности энергии поля е и потерь энергии поля Б ¿у:

де дБ ¿у

т - 1 _ .

Здесь I — средний по периоду волны вектор Пойнтинга,

IV I + , .2 I р2 I | т| т| т|

16пе _ |Е2| + |Н2| + V |Ут|„ ^|Рт| +

^ Вти2т

С „2„ (оГл М|р2|2 /г^0 / 2 , 2 , о__1т2

+-п2п2 3(1 - /г)|Е2|2 + (д2 + р2 + 2т^д|Е2| 8п \ тя

dDLV ттг ^ cn2n2 frRo , 2 2 , оIч

иг=<Ej>+W"E+«^т<р2+q2 - r-q|E2|)-

Потери энергии Dlv в итоге выделяются в виде тепла в среде и приводят к изменению физических свойств среды.

В настоящей работе использованы следующие параметры для плавленого кварца [1, 4]: ш = 0.64ше, Eg0 = 9 эВ, Ех = 0.6962, B2 = 0.4079, B3 = 0.8975, Ai = 0.0684, Л2 = 0.1162, A3 = 9.8962 мкм. Это даёт — = 27.539, -2 = 16.21, -3 = 0.19034 фс-1. При Л = 800 нм (— = 2.35 фс-1) имеем n = 1.45. Остальные параметры таковы: n2 = 2.48 ■ 10-16 см2/Вт, а = 6, WPIo = 3.7 ■ 1034 см-3с-1, po = 6.6 ■ 1022 см-3, rd = 32 фс, ts = 12 фс, Ttr = 150 фс, —тс = 3.

Для плавленого кварца fr = 0.18, однако в работе для сокращения трудоёмкости расчётов полагается fr = 0. По этой же причине в (9) пренебрегалось производной по времени d<pv)/dt. Тогда для потерь энергии имеем

nir = + wpi Ego(1 + |E2|/<4E2)). <17)

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

2. Нормировка

В исходных уравнениях перейдем к безразмерным переменным. Будем измерять время в t* = 1 фс, расстояние в x* = 1 мкм, частоты в фс-1, волновой вектор в мкм-1, скорость света c в мкм/фс, плотность в единицах начальной плотности атомов po, величины E, D, P, H в E*, определяемой формулой (16), скорость электронов v в выражениях (9), (11) в единицах v* = cE*/<4nex*po). Тогда (1) примет вид

-- i-D = pv + rot H - C1 <1 - p)|E2|a-1<1 + |E2|/4)E. <18)

c dt c

Уравнения (2)-(4), (6)-(8) и (12) в безразмерном виде выглядят так же, как и в размерном, а уравнения (5), (9) и (10) записываются как

Pn = C2 (<1 - fr )|E|2 + fr Rod) E, <19)

^ - i-pv = -C3PE - ^, <20)

dt Tc

I = C4<|E|2).<1 - p) + C- A <21)

При fr = 0 и d<pv)/dt = 0 уравнения (2), (20) и (17) выглядят следующим образом:

D = E + ^ Pm + C2|E|2E, <22)

m

C3TcE ioo\

v = "-7, <23)

i—Tc — 1

SDf = S (Т+S + Cd - Р)1Е21"(1 + ^/4)) • (24)

dt 8п V 1 + ш

Здесь C1 5 — соответствующие безразмерные постоянные.

3. Цилиндрически симметричный случай

Будем рассматривать задачи, в которых в цилиндрической системе координат зависимость от азимутального угла имеет вид

(Er, Ez, H^) = (Er, Ez, H^) cos (E^, H, H) = (E^, H, Я) sin <£•

Чтобы такая зависимость имела место, необходимо выполнение неочевидного равенства |Er|2 + |EEz|2 = |E^|2, которое предполагается в приближении уравнения Шрёдингера. В этом приближении Ez = 0, а уравнения для Er и E^ совпадают. Поэтому можно ожидать (и проверить в расчётах), что в случае уравнений Максвелла зависимость |E21 от угла будет слабой. Тогда для |E2| можно использовать выражение

|E2|

|E |2 + |E^|2 + |Ez |2 2

Использованное приближение будет верным, если |Ег|2 + |Ег|2 - |Е^|2 ^ |Е2|.

В случае обсуждаемой зависимости от угла величины с "шапочкой" описываются уравнениями (18)-(22), в которых операторы ротора и дивергенции в цилиндрической системе координат имеют вид

Ez

<9E^ dEr dEz 1 d(rE^) E

r dz ' dz dr ' r dr

+ — r

rot H =1 Я -

di di 1 d(rHv) Я

dz ' dz

dr ' r dr

div e = 1 dínEl + ^ + fi

r dr r dz

div H = -

1 d(rHr) Я , dHz

r dr

--^ +

r

dz

r

r

4. Начальные и граничные условия

В начальный момент времени £ _ - ¿о ^ - т_£ полагается, что в расчётной области отсутствуют электромагнитное поле (Е _ Н _ Рт _ \т _ 0) и ионизация (р _ 0). На оси г _ 0 имеем обычные для рассматриваемой моды граничные условия

Ег _ -Е^, Нг _ Н, Ег _ Н _ 0. (25)

На бесконечности наиболее адекватным является граничное условие Е _ 0, Н _ 0. В расчётах авторы вынуждены использовать ограниченную область. Условие Е _ 0,

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

пЕг — Ну = 0, пЕу + Нг = 0 при г = г0,

пЕу — Н = 0, гпЕх + = 0 при г = г0.

Тестовые расчёты показали, что такие граничные условия устраняют проблему отражённой волны, при этом их реализация значительно проще, чем, например, использование так называемых поглощающих граничных условий [9].

При г = 0 задается условие на инварианты Римана на входящей характеристике, соответствующее фокусируемому линзой пучку с линейной поляризацией

пЕг + = 2пЕ0, пЕу — Нг = — 2пЕ0, (26)

Ео = Етг ехр(—г2/т2 — ¿2/т| — г^0Г2/(2/)),

т = т/ (1 + & 2/4 )1/2, / = & + //&, г/ = /2. (27)

Здесь к0 = пш/с — волновой вектор в среде, т^ — длительность импульса, г/ — длина Рэлея, / — кривизна луча на расстоянии & от линейного фокуса, т/ = 1 мкм.

В пределе к0т ^ 1 (волна распространяется почти параллельно оси г) и при отсутствии волн, приходящих к границе г = 0 из расчётной области (т. е. при равенстве нулю инвариантов Римана на выходящей характеристике: пЕг — = 0, пЕу + Нг = 0), амплитуда волны Егпг связана с полной энергией импульса выражением (в размерных величинах)

Е

гпг

16£ш 4 1/2

пст2т^(п/2)1/2/

а условие (26) эквивалентно условию Ег = Е0, Еу = —Е0. При достаточно большом значении & этот предел действительно выполняется: в расчётах заведомо к0т > 100, а волны, отражённые от фокальных областей (расположенных на расстояниях г ~ &) вследствие нелинейности, приходят к границе г = 0 в моменты времени, когда величина Е0 на границе уже мала. Заметим, что условие Ег = Е0, Еу = —Е0 даёт практически те же результаты, что и (26), так как при достаточно большом & волне, отражённой от границы г = 0, требуется значительное время для возвращения в физически значимую фокальную область и амплитуда ее мала.

5. Конечно-разностная схема

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

Узлы сетки, совпадающие с границами расчётной области, назовем целыми, а сдвинутые относительно них на полшага — полуцелыми. Будем вычислять р в полуцелых узлах, а компоненты полей E и H в узлах Eri+i/2,fc, E^,fc, £¿¿^+1/2, Hri,fc+i/2, H^+i/2,fc+i/2, Hzi+i/2,fc. На такой сетке пространственные производные аппроксимируются центральными разностями со вторым порядком точности естественным образом. Например, z-компонента ротора от H, необходимая при вычислении z-компоненты электрического поля, имеет вид

, tjs 1 ri+i/2H^i+i/2,fc+i/2 — ri-i/2H^i-i/2,fc+i/2 Hr i,fc+i/2

(rot Hjj,k+i/2 —----.

Г» Ti+i/2 - fi-i/2 Г»

Поскольку в направлении z волна достаточно близка к монохроматической с волновым вектором k, то при вычислении производных в качестве шага по z использовалась величина hz — (2/k) sin(khz/2) « hz(1 + (khz)2/48). Такой выбор шага обеспечивает точное вычисление производных по z для параметров строго монохроматической волны, сохраняет второй порядок аппроксимации для остальных гармоник и даёт существенный эффект для точности получаемого решения.

Рассмотрим реализацию условий на оси в конечно-разностном случае. Будем вычислять E<£ и Hr при i — 0 сквозным расчётом, т. е. так же, как при i — 0. Для этого доопределим Hz в фиктивных точках Hzi=-i/2 — — Hzi=i/2. Для вычисления Hr на оси используем то обстоятельство, что z-компонента rot H на оси равна нулю (что требует симметрия задачи и, в частности, уравнение для Ez на оси). Это даёт

1 д (rHy) Hr — 0

г дг г

Интегрируя данное выражение с весом г от 0 до ri=i/2, в соответствии с (25) получим граничное условие Hri=o,fc+i/2 — H^=i,fc+i/2.

Описанная пространственная аппроксимация обеспечивает точное выполнение конечно-разностного аналога равенства rot VF — 0. Можно показать, что достигается второй порядок аппроксимации по пространству с учётом условий на оси.

Для аппроксимации (27) используем выражение

nErfc=M+i/2 + Er fc=2,i+i/2 + H — 2nE (t Г )

n-2--+ H^fc=i/2,i+i/2 — 2nE0(i, ri+i/2).

Условия на характеристиках на остальных границах будем аппроксимировать аналогично. Производные по времени d//dt аппроксимируем выражением (fn+i — /п)/т (т — шаг по времени), а члены уравнений вида iw/ — iw(/n+i + /n)/2.

Поскольку нас интересует процесс распространения волны, то ограничение на шаг по времени типа условия Куранта ст < h не является обременительным и позволяет использовать явную схему. Однако возможны такие параметры сетки, при которых ограничение вида wmT < 1 окажется более жёстким, чем условие Куранта. Поэтому уравнения (18), (22), (3), (4) аппроксимируются по времени следующим образом (F — несущественная в данном случае правая часть (18)):

ТЛп+i pjra T»n+i _L "Л™

D-— — iwD n+ D — cF™, (28)

т2

Dn+i — En+i + J] Pm+i + C2|En|2En+i, (29)

pn+1 _ pn pn+1 I pn Vra+1 + Vn

— гш---=---, (30)

T 2 2

Vn+l _ Vn Vn+1 + Vn /pn+1 I pn En+1 + E

m m-,,m'm ,.2 1 m ' m i-> i ■L-'

• . m 1 m _ m 1 -*- m ту __' 1 Coi ^

T — -2- = —1^m l -2---2- /' ()

Для решения системы (28)—(31) положим

pm+1 = + CPmEn+1, vm^1 = Cvm + CvmEn+1. (32)

Коэффициенты (рт, ^Ут, (ут находятся подстановкой (32) в (30) и (31). Затем, зная Юга+1 из (28) и подставляя (32) в (29), вычисляем Ега+1. После этого по формулам (30), (31) находим Рт+1, Ут+1.

Таким образом, схема содержит следующие этапы реализации:

1) вычисляются значения |ЕП|2 в полуцелых точках. В дальнейшем эти величины используются при вычислении нелинейных членов;

2) вычисляется значение рга+1;

3) вычисляются значения поля Юга+1, при этом величины |ЕП|2 и рга+1 интерполируются линейным образом из полуцелых узлов в узлы вычисления компонент Ю;

4) рассчитываются Р^^1, ^т+1, Ега+1 во внутренних точках;

5) задаются значения Ега+1 на границах;

6) по значениям Ега+1 вычисляется Ыга+1.

Представленная схема имеет второй порядок аппроксимации по пространственным переменным и первый — по времени. В линейном приближении при надлежащем задании начальных условий данная схема имела бы второй порядок аппроксимации по времени [16].

Приведенная конечно-разностная схема легко модифицируется при / = 0,

5 = 0.

6. Методические расчёты

Рассмотрим влияние различных параметров поставленной выше задачи на точность численного решения. Сравним результаты исследования для двух импульсов. В первом длительность импульса tl = 80 фс, энергия £1N = 2.5 мкДж, во втором tl = 150 фс, £in = 1 мкДж. Отметим, что такие условия типичны для экспериментов по лазерной модификации плавленого кварца [1-4, 6].

Приведем результаты методических расчётов для импульса tl = 80 фс.

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

2. Размер расчётной области по радиусу зависит от параметра d (27) (расстояние от границы до геометрического фокуса). Первоначальный радиальный размер лазерного импульса w пропорционален d. Расчеты показывают, что при d =120 достаточно ограничиться размером области по радиусу г0 = 60 мкм.

При фокусировке пучка электромагнитное поле существенно отлично от нуля на расстояниях от оси менее 10 мкм, а плотность плазмы — на расстояниях нескольких микрометров. Поэтому для временного интервала, когда лазерный импульс находится в зоне фокуса, была предпринята попытка существенно уменьшить размер расчётной области по г, сократив таким образом время решения задачи. Однако, как показывают расчёты, этого делать нельзя, поскольку при таком уменьшении г0 с течением времени решение грубо искажается, что, по-видимому, связано с наличием волн, сходящихся к оси с увеличивающейся амплитудой. Сильное уменьшение расчётной области исключает эти волны из рассмотрения. Незначительное же уменьшение расчётной области не имеет смысла, так как не приводит к существенному сокращению времени расчёта. По этой же причине увеличение шага сетки по г на периферии также оказывается не рациональным. На рис. 1 приведена зависимость ртах = тах р от времени при разных

т,х

шагах равномерной сетки и в случае использования неравномерной по г сетки с шагом Нг = 0.0707 мкм при г < 10 мкм, плавно увеличивающемся до 0.1 мкм к периферии. Видно, что кривая в случае неравномерной сетки достаточно далеко отклоняется от варианта с постоянным шагом Нг = 0.0707 мкм в сторону варианта с Нг = 0.1 мкм.

3. Шаги конечно-разностной сетки. Из рис. 1 наблюдается сходимость численного решения с уменьшением шага сетки. Для демонстрации сходимости выбрана плотность свободных электронов р, поскольку она наиболее чувствительна к шагу сетки. Это связано с тем, что скорость генерации свободных электронов пропорциональна |Е|12, а следовательно, небольшие изменения в величине электрического поля при изменении шага сетки приводят к значительным изменениям в распределении р. Из рис. 1, а также анализа других результатов расчётов можно сделать заключение, что точность, получаемая при Нг = Нг = 0.0707 мкм, вполне удовлетворительна.

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

Рис. 1. Зависимость pmax от времени в случаях равномерной по r сетки (hr = 0.1 мкм (1), 0.0707 мкм (2), 0.05 мкм (3)) и неравномерной по r сетки (4) при tl = 80 фс, eIN = 2.5 мкДж; 5 — тl = 150 фс, eiN = 1 мкДж

1.0

0 500 1000 t

Рис. 2. Зависимость E^ах от времени: 1, 2, 3 — tl = 80 фс, Ein = 2.5 мкДж, d = 160, 120, 80 мкм соответственно; 4 — tl = 150 фс, Ein = 1 мкДж, d = 120 мкм

4. Расчёты показали, что схема устойчива при т < const/(h-1 + h-1), где const ^ 0.5. Погрешности, связанные с конечным значением т, намного меньше погрешностей, связанных с конечным шагом сетки по пространству.

5. Зависимость от параметра d. На рис. 2 представлена зависимость Е^ах = max |E2|

r,z

от времени при разных значениях параметра d — расстояния от границы расчётной области до геометрического фокуса. На рисунке для наглядности сделан сдвиг на соответствующую разность по времени при разных d, необходимую пучку для достижения фокуса. Видно, что при d > 120 мкм различия в поведение Е^ах невелики. Исключение составляет Е^ах в начальные моменты своего роста. В отсутствие упомянутого сдвига для разных d эти моменты соответствуют временам t < 0, когда максимум |E2| достигается на границе z = 0. Для наших исследований подобные детали интереса не представляют. Заметим также, что при меньших d можно брать меньший размер области по радиусу, поскольку необходимо, чтобы величина r0 была существенно больше w пропорционально d.

6. Расчеты показали, что отклонения от цилиндрической симметрии действительно малы. Величина |Er|2 + |Ez|2 — |E^|2 не превышает 2-3% от |E2|.

При лазерном импульсе с Ein = 1 мкДж и tl = 150 фс результаты методических расчётов отличаются от рассмотренного случая Ein = 2.5 мкДж, tl = 80 фс только тем, что для достижения достаточной точности расчётов необходимо брать больший размер области по z (300 мкм) и можно ограничиться меньшими значениями d (80 мкм).

7. Физические результаты

Распределения |E21 и плотности для пучка с параметрами tl = 80 фс и Ein = 2.5 мкДж в различные моменты времени представлены на рис. 3 и 4. Согласно расчётам, картина распространения импульса в этом случае следующая.

Рис. 3. Распределение |E21 (слева направо) t = 350, 450, 500, 550, 600 фс. Максимальные значения |E2| при этом составляют 0.455, 0.882, 0.683, 0.42, 0.49 соответственно. Импульс с параметрами tl = 80 фс, Ein = 2.5 мкДж

Рис. 4. Распределение плотности р (слева направо) в моменты времени t = 410, 500, 550, 600 фс (ртах = 0.0025, 0.0081, 0.0059, 0.008 соответственно) и поверхность р при t = 600. Импульс с параметрами tl = 80 фс, Ein = 2.5 мкДж

На первом этапе импульс испытывает обычную фокусировку, которая благодаря эффекту Керра происходит на расстояниях от границы z = 0, несколько меньших геометрического фокуса. Фокусировка импульса приводит к образованию плазмы, поглощающей энергию электромагнитной волны. Появление плазмы вызывает уменьшение напряжённости поля в области, занимаемой плазмой, и вытеснение области с максимальным значением |E2| на периферию. В свою очередь это приводит к возникновению области высокой концентрации плазмы уже вне оси и способствует (наряду с силь-

2. МКМ МКМ

Рис. 5. Распределение энерговклада Dlv (Дж/см3) при tl = 80 фс, Ein = 2.5 мкДж (a) и TL = 150 фс, Ein = 1 мкДж (б)

ной зависимостью скорости ионизации от |E2|) образованию резкой границы занятой плазмой области с формированием больших градиентов р по радиусу. В некоторые моменты времени глобальный максимум |E2| достигается не на оси системы (см. рис. 3, t = 500). Распределение р также имеет большой локальный максимум вне оси (см. рис. 4, t = 550). Отметим высокую, превышающую в максимуме 1020 см-3, плотность плазмы. Характерные размеры, на которых происходят изменения различных величин, имеют порядок 1 мкм, в том числе по z, что нарушает условие применимости уравнения Шрёдингера. Отметим также наличие мощной рассеянной волны (см. рис. 3, t = 500 и 550), что не характерно для НУШ.

На временах, когда интенсивность лазерного импульса уже достаточно слаба, происходит отрыв волны от плазмы, плотность электронов испытывает простое экспоненциальное затухание со временем Ttr, импульс приобретает обычную форму и дефокуси-руется.

Область большого энерговклада Dlv (рис. 5) примерно соответствует области повышенной плотности электронов и имеет те же рассмотренные выше характерные особенности, хотя точного совпадения этих областей нет. Максимальное значение Dlv достигает нескольких тысяч Дж/см-3.

Таким образом, полученная зависимость максимальных по времени и пространству значений энерговклада, плотности свободных электронов и напряжённости электрического поля от Ein является не очень сильной (см. рис. 1, 2, 5). Наблюдается так называемый эффект клампинга [6]. Увеличение мощности импульса приводит к тому, что образование плазмы и потери энергии волны возникают раньше. При этом увеличивается пространственный размер области, в которой р, Dlv имеют большие значения.

Заключение

В работе предложена модель, описывающая распространение мощного фемтосекунд-ного лазерного импульса в нелинейной прозрачной среде с учётом генерации плазмы и дисперсионных свойств среды на основе уравнений Максвелла. Необходимость раз-

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

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

[1] Couairon A., Sudrie L., Franco M. et al. Filamintation and damage in fused silica induced by tightly focused femtosecond laser pulses // Phys. Rev. B. 2005. Vol. 71. Paper 125435 (11 p.).

[2] Couairon A., Mysyrowicz A. Femtosecond filamentation in transparent media // Phys. Rep. 2005. Vol. 441. P. 47-189.

[3] Winkler S.W., Burakov I.M., Stoian R. et al. Transient response of dielectric materials exposed to ultrashort laser radiation // Appl. Phys. A. 2006. Vol. 84. P. 413-422.

[4] Burakov I.M., Bulgakova N.M., Stoian R. et al. Spatial distribution of refractive index variations induced in bulk fused silica by single ultrashot and short laser pulses //J. Appl. Phys. 2007. Vol. 101. Paper 043506 (7 p.).

[5] Mermillod-Blondin A., Burakov I.M., Bulgakova N.M. et al. Flipping the sign of refractive index changes in ultrafast and temporally shaped laser-irradiated borosilicate crone optical glass at repetition rate // Phys. Rev. B. 2008. Vol. 77. Paper 104205 (8 p.).

[6] Булгакова Н.М., Стоян Р., Розенфельд А. Лазерно-индуцированная модификация прозрачных кристаллов и стекол // Квантовая электроника. 2010. Т. 40. С. 966-985.

[7] Turitsyn S.K., Mezentsev V.K., Dubov M. et al. Sub-critical regime of femtosecond inscription // Optics Express. 2007. Vol. 15, No. 22. P. 14750-14764.

[8] Ландау Л.Д., Лифшиц Е.М. Электродинамика сплошных сред. М.: Наука, 1982. 444 с.

[9] Petropoulos P.G. Reflectionless sponge layer for the numerical solution of Maxwell's equations in cylindrical and spherical coordinates // Appl. Numerical Math. 2000. Vol. 33. P. 517-524.

[10] Yuqing Jiao, Yaocheng Shi, Daoxin Dai, Sailing He. Accurate and efficient simulation for silicon-nanowire-based multimode interference couplers with a 3D finite-element mode-propagation analysis // J. Opt. Soc. Amer. B. 2010. Vol. 27, No. 9. P. 1813-1818.

[11] Odajima W., Tawa F., Shin-ya Hasegawa. Optical and thermal simulator for laser-assisted magnetic recording // Optical Rev. 2007. Vol. 14, No. 4. P. 180-185.

[12] Madsen N.K., Ziolkowski R.W. Numerical solution of Maxwell's equations in the time domain using irregular nonorthogonal grids // Wave Motion. 1988. Vol. 10. P. 583-596.

[13] Taflove A., Hagness S.C. Computational Electrodynamics: The Finite-difference Timedomain Method. 3rd edit. Artech House, Norwood, MA, 2005. 1038 p.

[14] Прокопьева Л.Ю., Федорук М.П., Лебедев А.С. Параллельный алгоритм метода конечных объёмов для решения трёхмерных уравнений Максвелла в нанокомпозитных средах // Вычисл. методы и программирование. 2009. Т. 10, № 1. С. 32-37.

[15] Popov K.I., McElcheran C., Briggs K. et al. Morfology of femtosecond laser modification of bulk dielectrics // Opt. Express. 2010. Vol. 19. P. 271-282.

[16] Федорук М.П., Березин Ю.А. Моделирование нестационарных плазменных процессов. Новосибирск: Наука, 1993. 356 с.

Поступила в 'редакцию 13 февраля 2012 г.