CC BY
0
Модель распространения световой пули в кристалле LiF
А.В. Кузнецов1
1 Иркутский филиал ФГБУН Института лазерной физики СО РАН, 664033, Россия, г. Иркутск, ул. Лермонтова, д. 130а
Аннотация
Построена модель распространения световой пули в кристалле LiF и формирования трека пули, состоящего из центров окраски. Модель воспроизводит наблюдавшиеся ранее продольные периодические колебания концентрации центров окраски в треке. Кроме того, модель впервые воспроизводит экспериментально наблюдаемые периодические колебания ширины трека. Относительный сдвиг между указанными колебаниями по результатам моделирования соответствует экспериментальным наблюдениям.
Ключевые слова: световая пуля, филаментация, центр окраски, LiF, FDTD.
Цитирование: Кузнецов, А.В. Модель распространения световой пули в кристалле LiF / А.В. Кузнецов // Компьютерная оптика. - 2024. - Т. 48, № 5. - С. 669-675. - DOI: 10.18287/2412-6179-CO-1419.
Citation: Kuznetsov AV. Model of light bullet propagation in LiF crystal. Computer Optics 2024; 48(5): 669-675. DOI: 10.18287/2412-6179-CO-1419.
Введение
Высокоинтенсивные фемтосекундные лазерные импульсы испытывают самофокусировку в прозрачных средах, формируя световые нити - филаменты. В твердых материалах филаменты могут оставлять треки с модифицированными оптическими свойствами. Этот эффект нашёл применение для формирования оптических элементов, таких как волноводы и дифракционные решетки, в объёме материалов. С научной точки зрения, формирование треков может использоваться для исследования взаимодействия высокоинтенсивного света и вещества.
При отрицательном параметре дисперсии групповой скорости д2к/ дю2 (здесь к - волновое число, ю -циклическая частота излучения) [1] реализуется особый случай филаментации, при котором световой импульс не только самофокусируется, но и сжимается в продольном направлении из-за совместного действия фазовой самомодуляции и дисперсии. Сжатый таким образом импульс называется световой пулей (СП) [2].
Удобным материалом для изучения формирования и распространения СП является кристалл фторида лития ЫБ, в котором филаменты оставляют хорошо видимые треки, состоящие из стабильных люминесцентных центров окраски (ЦО), таких как Б2 и Б3+ ЦО [3]. В этом смысле ОБ является нелинейным фоточувствительным материалом с высокой разрешающей способностью.
В 2015 году был проведен эксперимент с формированием СП в ОБ [4, 5]. Благодаря отмеченной фоточувствительности ОБ обнаружен эффект «дыхания» СП - периодической осцилляции пиковой напряженности электрического поля в СП при её движении. Экспериментально дыхание проявилось в периодической продольной осцилляции интенсивно-
сти люминесценции треков, наведенных СП. Дыхание СП получило универсальное объяснение дисперсионной природы, не связанное с какими-либо уникальными особенностями ОБ. Впоследствии дыхание СП наблюдалось по свечению плазменных каналов при прохождении СП в СаР2 [6], где не остаётся заметных треков. При исследовании люминесцентных треков СП в ОБ обнаружено, что ширина поперечного профиля трека (измеренная на половине высоты) испытывает продольные осцилляции практически с таким же периодом, как интенсивность люминесценции [5]. Причём, как видно на рис. 1 (рис. 4 в [5]), имеет место фазовый сдвиг между данными колебаниями. Его можно оценить по относительному положению максимумов и минимумов величиной в пределах 15 - 25 % от периода колебаний, который составляет в данном случае 31 - 32 мкм. Этот сдвиг наблюдался во всех исследованных треках СП в рамках работы [5]. Если бы фазовый сдвиг составлял около 50 % от периода колебаний (колебания в противофа-зе), его было бы легко объяснить тем, что уменьшение радиуса СП сопровождается одновременным увеличением напряжённости электрического поля. Компьютерное моделирование распространения СП в Ь^ показало именно такую картину — колебания радиуса СП и пиковой напряженности поля находятся в противофазе [7] (модель основана на уравнении однонаправленного распространения импульсного излучения). Таким образом, наблюдаемый эффект не поддается тривиальному объяснению и не был воспроизведен в компьютерном моделировании, судя по опубликованным работам. Чтобы приблизиться к пониманию эффекта, в данной работе предпринимается попытка моделирования распространения СП с сопутствующим наведением трека в Ь^ при помощи метода РБТБ, относительно свободного от аппроксимаций. Успешное воспроизведение обсуждаемого
эффекта проясняет, какой набор физических аспектов взаимодействия света и вещества необходимо учитывать для воспроизведения эффекта.
Готовых моделей на основе РБТБ-метода, учитывающих предполагаемый необходимым набор физических факторов (дисперсию, кубическую нелинейность, туннельную фотоионизацию, наведение трека) и располагающих явным описанием всех уравнений, необходимых для воспроизведения расчётов, в публикациях не найдено. В [8] и [9] обсуждаются РБТБ-модели с учетом нелинейной фотоионизации. Однако в них не учтена дисперсия и кубическая нелинейность; итоговые конечно-разностные уравнения не записаны. В [10] учтена дисперсия и различные виды нелинейности среды, но не учитывается ионизация. В
1 4-,
[11] учтены дисперсия, кубическая нелинейность, многофотонная фотоионизация, однако не записана итоговая система конечно-разностных уравнений. То же можно сказать о работе [12], где не приведены конечно-разностные уравнения для ионизации.
С учётом изложенного, в данной работе составляется и подробно описывается модель с необходимым набором физических факторов. Поскольку модель основана на общих принципах (не специализирована для воспроизведения одного эффекта), она может оказаться полезной для будущих исследований формирования и распространения СП в различных условиях, а также для приложений в области формирования оптических элементов в объёме прозрачных диэлектриков.
г 4 ^ о
3,5
3 -
2,5
3,5
- 3
2,5
20
40 мкм
60
н я 5
Рис. 1. [5] Внизу: изображение люминесцентной структуры, наведённой в LiF единичным 100-фс 3000-нм лазерным импульсом, распространявшимся слева направо; вверху: зависимости осевой интенсивности люминесценции (сплошная кривая 1) и диаметра структуры (пунктирная кривая 2) от продольной координаты
1. Физико-математическая модель филаментации 1.1. Уравнения Максвелла
Распространение светового импульса моделируется при помощи уравнений Максвелла
УхЕ = -ЫВ, Ы
с 2 ухв=ЫЕ3
Ы Во
(1) (2)
Здесь Е(г, ¿) - напряженность электрического поля, В(г, ^ - индукция магнитного поля, 3(г, {) - общая плотность тока, связанная с движением как связанных, так и свободных зарядов, о - скорость света, Во -электрическая постоянная.
Взаимодействие света со средой задается при помощи слагаемых, входящих в 3(г, ¿):
3 = 31 + 3 к + 3 а.
(3)
Здесь 3(г, ¿) - плотность тока, связанная с движением гармонических электронных осцилляторов в рамках модели рефракции Лоренца; 3(г, ^ - плотность тока, вносящая в модель мгновенный эффект Керра; 3а(г, {) — плотность тока, соответствующая нелинейному поглощению света.
1.2. Генерация электронных возбуждений
Одновременно с распространением импульса происходит генерация электронных возбуждений (ЭВ) -экситонов и электронно-дырочных пар. Их концентрация Ы(г, ¿) увеличивается в каждой точке пространства г под действием электрического поля согласно уравнению
дЫ
= Ж (|Е|)
1 -
N
Ыт
(4)
Здесь функция Ж(|Е|) описывает зависимость скорости роста концентрации ЭВ N от модуля напряженности электрического поля. Множитель 1 - N/ Ытах введен для ограничения роста концентрации ЭВ максимальной величиной Ктах = 6,117х1028 м - 3, равной концентрации ионных пар Ы + f - в кристалле ЫБ. Данная величина найдена как отношение плотности ЫБ р = 2635 кг / м3 [13] к сумме атомных масс Ы и Б:
Ытах = р / (ти + тР).
Параметр Келдыша у = ©Vти / еЕ0 [14] имеет величину 0,72 для монохроматического излучения частоты © = 6,08х1014 с-1, соответствующей длине волны X = 3,1 мкм, ширине запрещенной зоны и= 1,18x10 -18 Дж = 13,6 эВ [13] и амплитуде колеба-
ний напряженности электрического поля Е0 = 7,44*109 В / м, соответствующей интенсивности I = спеЕо / 2 = 1013 Вт / см2 при показателе преломления п = 1,365. Неравенство у < 1 означает, что туннельная ионизация преобладает над многофотонной. Поэтому для описания скорости генерации ЭВ в единице объёма под действием электрического поля напряженности |Е| в данной работе используется модель туннельной ионизации. В [14] приведена формула (40) для скорости туннельной ионизации ^(Е0), усредненной по периоду гармонического колебания электрического поля амплитуды Е0. Из данной формулы, с учётом нижеследующего комментария [14], выводится выражение для мгновенной скорости ионизации:
Ш (|Е|) = а.-
2 е | Е
и V 3п Й
2 ( ехр
V
пл/ти
~2 ей | Е |
Л
(5)
где а = 2,1 - множитель, величина которого подобрана так, чтобы приблизительно выполнялось равенство
1 г
М!(Е0)= — I Ш(ЕцБШю/)d(ю/).
тт •'0
(6)
1.3. Нелинейное поглощение
Нелинейное поглощение света и соответствующее рассеяние описывается плотностью тока
^ = и-
Е дЫ
| Е2 "дГ'
(7)
Данное уравнение выводится из равенства мощностей JaE = идЫ / д/ с учётом сонаправленности векторов ^ и Е.
1.4. Линейный показатель преломления
Линейный показатель преломления вводится в модель при помощи модели электронных осцилляторов Лоренца. Смещение осцилляторов от равновесных положений приводит к поляризации диэлектрика Р/ и соответствующей плотности тока J/ = дР/ / д/. Смещение Дг единичного осциллятора заряда е и массы т под действием электрического поля Е описывается уравнением
д 2Дг е -+ ю2 Дг = — Е,
д/2 т
(8)
где ю — собственная частота осциллятора.
При концентрации осцилляторов р поляризация диэлектрика равна Р/ = реДг. Из (8) следует уравнение для Р/:
дР+ю 2 Р/ = Р—
д/2 т
Е.
(9)
Вводя обозначение а2 = ре2 / те0 и индекс I для обозначения видов электронных осцилляторов с различными параметрами ю и а, перепишем (9) в виде системы уравнений
дJ„
+ ю2 РИ = а^Е,
д/
дР; Jl■ =—,
1 д/
J' = Б 4.
(10)
(11) (12)
Численные значения частот ю, для (10) найдены из данных, приведенных в работе [15]. Экспериментальная зависимость линейного показателя преломления п фторида лития от длины волны света в вакууме X аппроксимируется формулой Селлмейера
2(^) = 1+Б
а, X2 X2 — X2
(13)
где индекс , принимает значения 1 и 2, коэффициенты
а1 = 0,92549, а2 = 6,96747, резонансные длины волн в вакууме X1 =0,07376 мкм и X2 = 32,79 мкм. Учитывая соотношение
ю = -
2пс
Т"'
(14)
получаем значения резонансных частот со I для уравнения (9): ой! =2,55376 х1016 с -1, ¿52 = 5,74 4 5 9 х1013 с -1.
Согласно модели электронных осцилляторов Лоренца, зависимость показателя преломления от частоты описывается уравнением [1]
2(ю) = 1+Б
а,
, со2 — ю2
(15)
Приравнивая соответствующие слагаемые из правых частей (13) и (15) и учитывая (14), находим значения параметров а, для (10):
а2 = а1 ю2.
(16)
1.5. Нелинейный показатель преломления
Нелинейный показатель преломления вводится в модель при помощи плотности тока
Jк (г, /) = ,
д/
(17)
где поляризация Рк пропорциональна кубу мгновенной напряженности электрического поля:
Рк = 80Х(3)|Е|2 Е.
(18)
Здесь х(3) = 4,7^10 -23 м2/ В2 - скалярная диэлектрическая восприимчивость третьего порядка. Данное значение получено переводом экспериментального
п
п
значения x(3) = 0,34x10 -14см3/ эрг [16] из СГС в СИ при помощи соотношения xS? = 4я(3;
;104)_2 xg [1].
1.6. Метод конечных разностей во временной области
Для решения уравнений Максвелла в этой работе применяется метод конечных разностей во временной области (FDTD) [17, 18]. Из [18] позаимствован метод дополнительных дифференциальных уравнений (auxiliary differential equation method (ADE), раздел 9.4), применяемый в рамках FDTD для описания линейной (не зависящей от амплитуды колебаний поля) дисперсии в рамках FDTD-метода (раздел 9.4.2, уравнение (9.50)) и кубической нелинейности (9.6.4, уравнение (9.95)). В [18] также описан альтернативный метод для описания дисперсии - метод кусочно-линейной рекурсивной свертки (piecewise-linear recursive-convolution method (PLRC), раздел 9.3). Однако авторы [18] отмечают, что метод ADE удобнее, чем PLRC, при том же втором порядке точности.
Метод ADE позаимствован с определенными изменениями по сравнению с [18]. Уравнение (9.50) [18] выражает плотность тока в «промежуточный» момент времени как среднее арифметическое плотностей тока в соседних «целых» моментах. Чтобы избежать необходимости этого усреднения при описании дисперсии и аналогичного усреднения при описании нелинейной фотоионизации, дискретная пространственно-временная сетка в данной работе модифицирована по сравнению со стандартным FDTD-методом, как показано ниже.
Для снижения объема вычислений рассматривается двумерная задача c электромагнитной волной TMz-поляризации. Волна распространяется вдоль оси х, имея поляризацию вдоль перпендикулярной оси z. Поля волны постоянны вдоль оси z (пучок бесконечно широкий по z). Соответственно, все производные по z равны нулю. Тогда у векторов E, P, J отличны от нуля только z-компоненты, а у векторов B - х- и y-компоненты. Поперечное измерение y обеспечивает возможность самофокусировки и дифракции.
Запишем уравнения Максвелла (1) и (2) в скалярной форме для декартовой системы координат (х, y, z)
дЕ.
{zxy}
дЕ,
{yzx} _
дБ{
{xyz}
д{ yzx} d{zxy}
c2 | 9Bizxy}
дБ,
,{ yzx}
д{ yzx} д^>
дt
= дE{xyz} + J{xyz}
дt
(19)
(20)
Здесь выражение вида (хух) означает циклическую перестановку х^у^х^...
С учётом ТМ7-поляризации волны, уравнения (19, 20) упрощаются до
дБх дEz _ дБy
дЕг
дy дt дx дt
(21)
„2, дБу
дБ^ j= SEl + J дx дy I дt е0
(22)
Индекс z у векторных компонент Ez, Pz и Jz далее для краткости опускается.
Как и в традиционной FDTD-схеме, электрическое и магнитное поля в этой работе определены в пространственных узлах прямоугольной дискретной сетки со взаимным чередованием. Электрическое поле Е определено в точках xx = XAx, yY = YAy, где X и Y - целочисленные индексы, а Ax и Ay - пространственные шаги дискретной сетки. В этих же точках определены величины J, P и N. Магнитное поле определено в «промежуточных» точках: компонента Б;с - в точках xx = XAx, yY = (Y+1 / 2)Ay, а компонента Бу - в точках xx = (X +1 / 2) Ax, yY = YAy.
Искомые будущие пространственные распределения полей в FDTD-методе определяются методом перешагивания (leapfrog method), чередуясь во времени:
E(t = TAt) ^ Б(/ = (T +1/ 2) At) ^ E(t = (T +1) At)...
Здесь T - целое число, At - временной шаг. Внесём в схему модификацию, чтобы избежать упомянутой выше необходимости усреднения плотности тока между двумя моментами времени, как в уравнении (9.50) [18]. Добавим «параллельную» последовательность временных шагов на той же пространственной сетке
Б(/ = TAt) ^ E(t = (T +1/2) At) ^ Б(/ = (T +1) At)...
Теперь каждое из полей (электрическое и магнитное) определено как в «основные» моменты времени (кратные At), так и в «промежуточные» моменты. Поскольку эти моменты времени принципиально не различаются (содержат оба поля), устраним формальное различие между ними. Для этого введём вдвое меньший временной шаг: Ax = At /2. Тогда в формулах остаются только целые множители при временном шаге:
E(t = TAx) ^ Б(/ = (T +1) Ax) ^ E(t = (T + 2) Ax)...
Б(/ = TAx) ^ E(t = (T +1) Ax) ^ Б(/ = (T + 2)Ax)...
Описанная модификация сохраняет второй порядок точности FDTD-схемы.
Конечно-разностные аппроксимации уравнений (21, 22) представим следующим образом (верхние индексы обозначают расположение соответствующих величин в пространственно-временной сетке):
EX, Y,T _ EX
AT
EXYT _ E
Бх ,Y + 1/2, t + 1 _ Б
1X, Y + 1/2,t _1
2Ax
'X + 1, y,t БХ + 1/2,y,t + 1 _ Бх + 1/2,y,t_'
Ax
2Ax
+ 1/2,y,t _ БХ_1/2,y,t Бх,Y+1/2,t _ БХY_1/2,t
Ax
Ex ,Y ,t+1 _EX ,Y ,T _1 JX ,Y ,T
-+-.
Ay
(23)
(24)
(25)
2Ax
и
0
В уравнениях (23 - 25) искомыми величинами являются «будущие» значения компонент полей с временными индексами 7+1. Искомые компоненты магнитного
поля Б,
x ,y+1/2,t+1
и Б
X + 1/2, y ,t + 1
выражаются из (23, 24).
Искомое электрическое поле Ех,г,т+1 определяется из (25). Для решения этого уравнения необходимо предварительно найти входящую в него плотность тока
3 х т = 3 * -т + 3* т + 3* -т. (26)
Слагаемое 3х,т,т, согласно (4) и (7), равно
TTW( EXYT) L NX
EX
Nm
(27)
Фигурирующая здесь концентрация электронно-дырочных пар Ыхг,т выражается из конечно-разностной аппроксимации уравнения (4):
NX,Y,T _ NX,Y,T_2 I N
■ = W (E) I 1 _
2Дт
Nm
(28)
Слагаемое JX,Y,T определяется согласно уравнениям (17 - 18)
JX, Y,T J k
S0%
(3)
2Дх
|(EX Y T+1 )3 _(EXYT_ )3 ). (29)
Данное выражение содержит куб искомой напряженности электрического поля EX,Y,T+l, вследствие чего уравнение (25) оказывается кубическим.
Слагаемое Jf,Y,T состоит из i-х компонент, определяемых посредством решения системы уравнений (10, 11) методом с перешагиванием (leapfrog method). Конечно-разностные аппроксимации уравнений (10, 11) запишем следующим образом:
JX , y ,t _ JX, Y,T_ 2
2Дт
+ ю 2 P,XJ,T _' = a2 e0 EX
P,X
_1 _ dx, y ,t_ 3
1 /• /•
TX, y,t_2 _ _2_
Jh
2Дх
(30)
(31)
Из (31) выражается поляризация Р/7,т—1 и подставляется в (30). Далее из (30) определяется плотность тока 3Х/',т для ,-го типа электронного осциллятора. Наконец, вычисляется сумма:
3 х",т = БТ". (32)
Компоненты плотности тока, определенные согласно (27), (29), (32), суммируются для получения общей плотности тока 3х-г-т, которая, в свою очередь, подставляется в (25). Получается кубическое уравнение относительно искомой напряжённости электрического поля Ех,7,т+1
p3 (EXY-T+' j + pjEX•r•T+1 + p0 =0
(33)
с коэффициентами р3, р\ и р0 (явные выражения для коэффициентов не показаны для краткости). Это
уравнение решается итерационным методом Ньютона, причём в качестве начального приближения используется величина -р0 / р1, соответствующая решению (33) при р3 = 0.
Условие устойчивости описанной схемы - достаточно малый шаг по времени, подбираемый эмпирически.
1.7. Расчётная область, генерация импульса
Расчетная область имеет прямоугольную форму. Слева она ограничена гранью х = 0, с которой испускается линейно поляризованный импульс с плоским волновым фронтом при помощи задания граничного условия
Ez (x = 0, y, t) = A(y, t )sin(rot) с амплитудой гауссовой формы
A( y, t ) = A„exp
I -2 12 ^
d2y
d2
(34)
(35)
Здесь A0=25 В/м - пиковая амплитуда, dy=20 мкм и dt=40 фс - характеристики ширины и длительности импульса, ю - частота, соответствующая длине волны в вакууме 3,1 мкм. Величина амплитуды A0 здесь подобрана так, чтобы немного превышать минимально необходимую для самофокусировки. Начальные ширина и длительность импульса выбраны меньше, чем обычно реализуются в экспериментах для экономии вычислительных ресурсов, но достаточно большими, чтобы размеры импульса были значительно больше длины волны.
На правой грани x=xmax задано простое условие Ez = 0. Хотя это условие является отражающим, фактически отражения не происходит, так как расчёт останавливается прежде, чем импульс достиг бы этой грани. Величина xmax = 1,5 мм выбрана так, чтобы трек СП помещался в расчётной области.
Поскольку испускаемый световой импульс зеркально-симметричен относительно оси x, для оптимизации расчёты ведутся только с одной стороны от этой оси - при y > 0. Соответственно, на нижней грани расчётной области y = 0 задано симметричное граничное условие.
Чтобы ослабить отражение расходящегося из-за дифракции излучения обратно в расчётную область, желательно задать поглощающее граничное условие на верхней грани y = ymax. В этой работе применено условие типа «perfectly matched layer» (PML). Выбрана реализация PML, описанная в параграфе 3.2 монографии [19]. Потребовались тривиальные изменения её для адаптации к численной схеме в данной работе. Выбор обусловлен достаточно подробным описанием в указанном источнике.
Достаточно большая величина ymax = 10dy позволяет существенно ослабить краевые эффекты.
В начальный момент времени поля поляризация, плотности тока и концентрация электронных возбуждений равны нулю.
2. Результат моделирования
В результате моделирования получена временная эволюция электрического и магнитного полей импульса, а также пространственное распределение концентрации ЭВ Щх, у) в треке СП. Для простоты считается, что концентрация ЦО равна Ы(х, у), хотя в действительности лишь некоторая доля ЭВ трансформируется в стабильные ЦО.
На рис. 2 показано тоновое изображение распределения концентрации ЦО в треке Ы(х, у) и концентрация на оси трека, отнесенная к пиковой концентрации N (х) = N (х, у=0) / тах( N (х, у)).
2 ш ■
-2
Рис. 2. Тоновое изображение расчётного распределения концентрации ЦО в треке СП (вверху) и относительная концентрация ЦО на оси трека Щ(х) (внизу)
Видны колебания расчётной концентрации ЦО вдоль оси х, подобные наблюдавшимся в эксперименте. Период колебаний составляет около 30 мкм, что соответствует эксперименту [5].
На рис. 3 показан наиболее интенсивный фрагмент трека и соответствующие колебания концентрации ЦО на оси трека Щх) и ширины 5(х) поперечного профиля трека, измеренной на половине высоты. Колебания имеют практически одинаковый период, что соответствует эксперименту (рис. 1). Относительный сдвиг между колебаниями (точнее говоря, между их основными гармониками) составляет около 25 % от периода колебаний при том же направлении (колебания ширины смещены «вправо»), что приблизительно соответствует эксперименту (рис. 1).
При варьировании некоторых параметров в реалистичных пределах (ширины запрещённой зоны ЫБ, длины волны излучения, энергии импульса, начальных размеров импульса) относительный сдвиг колебаний изменялся на величину порядка единиц процентов от периода колебаний. Целенаправленного исследования в этом направлении не проводилось. Не было попыток найти допустимые пределы изменения параметров, за которыми результат моделирования оказывался бы явно нереалистичным. Также не производился подбор оптимальных параметров (подгонка) для приближения результата моделирования к эксперименту.
Таким образом, модель воспроизводит наблюдаемые колебания осевой плотности и ширины трека СП, а также сдвиг фазы между этими колебаниями.
Установлено, что колебания ширины трека СП в процессе его наведения совпадают по фазе с колеба-
ниями ширины центральной части СП, где напряженность электрического поля достаточно велика для заметной генерации центров окраски. Причём колебания ширины этой части СП смещены относительно колебаний пиковой напряженности поля в СП (на качественном уровне это пока не объяснено). Такое «поведение» поля «записывается» в треке в виде колебаний его ширины и интенсивности.
2 3
I о
^ — 3
1,75 1,0
1,50 о,5
Рис. 3. Расчётное распределение концентрации ЦО в треке СП (вверху); осцилляции концентрации ЦО на оси трека Ñ(x) и ширины трека 8(x) (внизу)
Заключение
Построена модель распространения фемтосекунд-ного лазерного импульса среднего ИК-диапазона с формированием из него световой пули в условиях аномальной дисперсии групповой скорости. Модель апробирована на примере распространения световой пули в кристалле LiF. При этом модель воспроизводит продольные колебания интенсивности люминесцентного трека из центров окраски, наведённого световой пулей. Впервые воспроизведены колебания ширины поперечного профиля трека, а также ранее наблюдавшийся в экспериментах сдвиг между продольными колебаниями ширины поперечного пофиля трека и интенсивности люминесценции.
Отметим преимущества данной модели по сравнению с моделями, обычно применяемыми для моделирования филаментации и световых пуль. В отличие от моделей на основе нелинейного уравнения Шрё-дингера или на основе уравнения однонаправленного распространения импульсного излучения, описанная модель на основе FDTD-метода не накладывает особых ограничений на относительные направления распространения излучения, на наличие рассеяния света средой и отражения. Соответственно, становится возможным моделирование распространения световых волн произвольной конфигурации в неоднородных средах. Таким образом, модель более приспособлена к техническим применениям.
Благодарности
Работа поддержана научным проектом 0243-20210004 в рамках плана фундаментальных исследований Российской академии наук на период до 2025 года.
References
[1] Boyd RW. Nonlinear optics. 4th ed. London: Academic Press; 2020. ISBN: 978-0-12-811002-7.
[2] Silberberg Y. Collapse of optical pulses. Opt Lett 1990; 15(22): 1282-1284. DOI: 10.1364/0L.15.001282.
[3] Martynovich EF, Dresvyansky VP, Rakevich AL, Lazareva NL, MA Arsentieva MA, Tyutrin AA, Bukhtsoozh O, Enkhbat S, Kostryukov PV, Perminov BE, Konyashchenko AV. Creating of luminescent defects in crystalline media by a scanning laser beam. Appl Phys Lett 2019; 114(12): 121901. DOI: 10.1063/1.5087688.
[4] Chekalin SV, Kompanets VO, Kuznetsov AV, Dormidonov AE, Kandidov VP. Regular 'breathing' of a near-single-cycle light bullet in mid-IR filament. Laser Phys Lett 2016; 13(6): 065401. DOI: 10.1088/16122011/13/6/065401.
[5] Kuznetsov AV, Kompanets VO, Dormidonov AE, Chekalin SV, Shlenov SA, Kandidov VP. Periodic colour-centre structure formed under filamentation of mid-IR femtosecond laser radiation in a LiF crystal. Quantum Electron 2016; 46(4): 379-386. DOI: 10.1070/QEL16038.
[6] Zaloznaya ED, Dormidonov AE, Kompanets VO, Chekalin SV, Kandidov VP. Material dispersion effect on the oscillations of a single-cycle wave packet. Optics Spectrosc 2022; 130(12): 1596. DOI: 10.21883/EOS.2022.12.55248.3933-22.
[7] Zaloznaya E, Kompanets V, Savvin A, Dormidonov A, Chekalin S, Kandidov V. Carrier-envelope phase effect on light bullet dynamics. Laser Phys Lett 2022; 19(7): 075402. DOI: 10.1088/1612-202X/ac7134.
[8] Zhang S, Tripepi M, AlShafey A, Talisa N, Nguyen HT, Reagan BA, Sistrunk E, Gibson DJ, Alessi DA, Chowdhury EA. Femtosecond damage experiments and modeling of broadband mid-infrared dielectric diffraction gratings. Opt Express 2021; 29(24): 39983-39999. DOI: 10.1364/OE.439895.
[9] Mahdy A. Numerical simulation for the efficiency of the produced terahertz radiation by two femtosecond laser pulses: Above-threshold-ionization. Journal of Applied
Mathematics and Physics 2023; 11(10): 2997-3008. DOI: 10.4236/jamp.2023.1110198.
[10] Tyrrell JCA, Kinsler P, New GHC. Pseudospectral spatial-domain: a new method for nonlinear pulse propagation in the few-cycle regime with arbitrary dispersion. J Mod Opt 2005; 52(7): 973-986. DOI: 10.1080/09500340512331334086.
[11] Bourgeade A, Mezel C, Saut O. Modeling the early ionization of dielectrics by ultrashort laser pulses. J Sci Comput 2010; 44: 170-190. DOI: 10.1007/s10915-010-9375-0.
[12] Schmitz H, Mezentsev V. Full-vectorial modeling of femtosecond pulses for laser inscription of photonic structures. J Opt Soc Am B 2012; 29(6): 1208-1217. DOI: 10.1364/JOSAB.29.001208.
[13] Weber MJ. Handbook of optical materials. CRC Press; 2003. ISBN: 0-8493-3512-4.
[14] Keldysh LV. Ionization in the field of a strong electromagnetic wave. Sov Phys JETP 1964; 20(5): 13071314.
[15] Li HH. Refractive index of alkali halides and its wavelength and temperature derivatives. J Phys Chem Ref Data 1976; 5(2): 329-528. DOI: 10.1063/1.555536.
[16] Levenson M. Feasibility of measuring the nonlinear index of refraction by third-order frequency mixing. IEEE J Quantum Electron 1974; 10(2): 110-115. DOI: 10.1109/jqe.1974.1145780.
[17] Yee K. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Trans. Antennas Propag 1966; 14(3): 302307. DOI: 10.1109/TAP.1966.1138693.
[18] Taflove A, Hagness SC. Computational electrodynamics: The finite-difference time-domain method. 3rd ed. Artech House Publishers; 2005. ISBN: 978-1-58053-832-9.
[19] Sullivan DM. Electromagnetic simulation using the FDTD method. The Institute of Electrical and Electronics Engineers; 2013. ISBN: 9781118459393.
Сведения об авторе
Кузнецов Андрей Викторович, 1982 года рождения, кандидат физико-математических наук по специальности «Лазерная физика», старший научный сотрудник Иркутского филиала ФГБУН Института лазерной физики СО РАН. Научные интересы: фемтосекундная лазерная филаментация (теория и эксперимент), центры окраски в кристаллах. E-mail: a.v.kuznetsov@bk.ru
ГРНТИ: 29.33.43 (распространение лазерного излучения).
Поступила в редакцию 30 августа 2023 г. Окончательный вариант - 24 декабря 2023 г.
Model of light bullet propagation in LiF crystal
A.V. Kuznetsov1 1 Irkutsk Branch of Institute of Laser Physics SB RAS, 664033, Irkutsk, Russia, Lermontov str. 130а
Abstract
A numerical model is developed for simulation of the propagation of a light bullet in LiF crystal and the formation of a colored track. The model reproduces the previously observed longitudinal oscillations in the concentration of color centers in the track and oscillations in the track width. For the first time, the phase shift between these oscillations is numerically reproduced.
Keywords: light bullet, filamentation, color center, LiF, FDTD.
Citation: Kuznetsov AV. Model of light bullet propagation in LiF crystal. Computer Optics 2024; 48(5): 669-675. DOI: 10.18287/2412-6179-CO-1419.
Acknowledgements: This work was financially supported by project № II.10.1.6 of the Siberian Branch of the Russian Academy of Sciences, "Mechanisms of extreme nondestructive interaction of solid dielectrics with intense laser radiation".
Author's information
Andrey Viktorovich Kuznetsov, (b. 1982) PhD of Physical and Mathematical Sciences in the field of Laser Physics. Currently he works as Senior Researcher at Irkutsk Branch of Institute of Laser Physics SB RAS. Research interests: femtosecond laser filamentation (both theory and experiment), color centers in crystals. E-mail: a.v.kuznetsov@bk.ru
Received August 30, 2023. The final version - December 24, 2023.
