Моделирование предвестников Зоммерфельда и Бриллюэна в среде с частотной дисперсией на основе разностного решения волнового уравнения

МОДЕЛИРОВАНИЕ ПРЕДВЕСТНИКОВ ЗОММЕРФЕЛЬДА И БРИЛЛЮЭНА В СРЕДЕ С ЧАСТОТНОЙ ДИСПЕРСИЕЙ НА ОСНОВЕ РАЗНОСТНОГО РЕШЕНИЯ ВОЛНОВОГО УРАВНЕНИЯ

Козлова Е.С., Котляр В.В.

Институт систем обработки изображений РАН

Аннотация

Волновое уравнение размерности (2+1) для ТЕ-поляризации с учётом частотной дисперсии материала записано в форме, в которой вторая производная по времени внесена под знак интеграла свёртки, описывающего электрическую индукцию. Предложен алгоритм решения такого уравнения с помощью явной разностной схемы. Решение волнового уравнения, описывающее распространение ультракороткого импульса длительностью 3,36 фс (в спектральном интервале 274-806 нм) внутри планарного волновода из кварцевого стекла на центральной длине волны 532 нм, отличается в среднем от решения этого уравнения без учёта дисперсии на 3% и от решения уравнений Максвелла, полученного FDTD-методом в программе FullWAVE, на 6%. Детальный анализ полученного решения для фемтосекундного импульса с резким начальным фронтом позволил обнаружить предвестники, которые приходят в точку наблюдения раньше основного импульса и по интенсивности в 100 раз меньше его. Причём рассчитанное время задержки предвестника Зоммерфельда отличается от теоретического на 20%, а предвестника Бриллюэна - на 3%.

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

Введение

В настоящее время большое количество экспериментальных и теоретических работ посвящено исследованию ультракоротких импульсов, длина которых соизмерима с длиной волны [1, 2]. Ультракороткие импульсы широко применяются в волоконной оптике [3], фемтосекундной и электронной микроскопии [4, 5], прецизионной микрообработке материалов [6, 7] и медицине [8]. Фемтосекундные импульсы применяют для генерации ультрафиолетового света [9] и в качестве «оптических пуль» [10].

Одним из фундаментальных свойств коротких импульсов, распространяющихся в среде с частотной дисперсией, является то, что высокочастотная составляющая импульса распространяется в среде со скоростью света в вакууме. Это явление было предсказано Зоммерфельдом и Бриллюэном в 1914 году. Данный феномен называется оптическим предвестником [11, 12]. Зоммерфельд и Бриллюэн асимптотическим методом седловой точки рассмотрели распространение синусоидального импульса, умноженного на функцию Хевисайда (резкий передний фронт импульса), в среде с одним резонансом Ло-рентца. Они нашли, что раньше всех в точку наблюдения приходит высокочастотный предвестник (предвестник Зоммерфельда - ПЗ), который распространяется со скоростью света в вакууме, следом за ним приходит другой предвестник (предвестник Бриллюэна - ПБ) с низкочастотной составляющей сигнала, который распространяется с фазовой скоростью сигнала, а за ним уже приходит основной сигнал с групповой скоростью. В [13] измерили скорость распространения в дисперсной среде переднего фронта короткого импульса, которая отличалась от скорости с света в вакууме всего на 0,01%. Моделирование распространения предвестников, как правило, проводят с помощью FDTD-метода, обобщённого на случай учёта дисперсии [14, 15]. В [14] рассматривался Ш FDTD-метод и моделировались предвестники в однополюсной

(один резонанс) дисперсионной среде Лорентца. В [15] моделировались предвестники в промежуточном спектральном режиме: несущая частота импульса не была равна частоте резонанса поглощения Лорентца и не была много меньше его, как в [14]. В [16, 17] предвестники исследуются аналитически также для однополюсной среды Лорентца. Показано, что амплитуда ПЗ описывается функцией Бесселя первого порядка, а ПБ - функцией Эй-ри. Заметим, что в работах [14-17] исследовались предвестники в модельной Ш дисперсной среде с одним резонансом Лорентца.

В данной работе моделируется распространение ПЗ и ПБ в 2D среде типа кварца с тремя резонансами в модели Селлмейера. Частотный режим использовался промежуточный, как в [15], но в отличие от [15] несущая частота была дальше от частоты резонанса. Обычно амплитуда предвестников очень мала и составляет 10-4 от максимальной амплитуды начального импульса. Поэтому для надёжного обнаружения предвестников (чтобы их амплитуда была больше амплитуды ошибки вычислений) мы подобрали в модели дисперсии Селлмей-ера резонансные частоты и амплитуды резонансов таким образом, что интенсивность ПЗ составляла 10-2 от интенсивности основного импульса и интенсивность ПБ соответственно 3 • 10-2.

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

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

ни диэлектрической проницаемости. Также существует модификация этого метода, в основу которой положено предположение о том, что функция диэлектрической проницаемости является кусочно-линейной на временном шаге [21]. Это позволяет повысить точность метода.

Ко второй группе относятся алгоритмы [22-28], использующие вспомогательное уравнение, связывающее плотность электрического поля и его напряжённость. Основным их преимуществом является высокая точность. Однако в дополнение к тем ресурсам, что использует стандартный FDTD-метод, требуется дополнительный объём памяти для хранения промежуточных переменных.

К третьей группе относятся алгоритм [29], в котором используется Z-преобразование уравнений Максвелла. Недостатком данного метода, как и в случае использования дополнительного дифференциального уравнения [26], является проблема ресурсоёмкости, высокой вычислительной сложности и временных затрат на получение решения.

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

т г к

временном слое только отсчёты напряжённости поля Е и отсчёт поля на предыдущем слое Е к-1. В FDTD-методе на каждом временном слое в случае ТЕ-поляризации хранятся отсчёты 3 проекций поля: Еу, Нх , Н2 . Проведено сравнение результатов численного моделирования, полученных с помощью разработанной схемы и программы FuИWAVE (RSoft), реализующей дисперсионный FDTD-метод. Моделирование распространения фемтосекундного импульса с резким передним фронтом в планарном волноводе с учётом дисперсии показало наличие у этого импульса ПЗ и ПБ, разделённых в пространстве.

1. Разностное решение волнового уравнения с учётом дисперсии

Распространение света в двумерной линейной изотропной диспергирующей среде для напряжённости электрического поля описывается волновым уравнением [11]:

а хА - Л =о,

х2 у л2

С

дг2

д2 д2

где Ах = —- +--- - оператор Лапласа, ЕУ - проекция

дх д2

вектора напряжённости электрического поля на ось у (ТЕ-поляризация); х и 2 - пространственные координаты; г - время; с - скорость света в вакууме; Dу -

проекция вектора электрической индукции на ось у. В случае немагнитной среды электрическая индукция будет представлена интегралом свёртки Дюамеля [11]:

D (х,2,г) = |е(х,2,?)Еу (х, 2,г -?)й?.

Воспользуемся принципом причинности (соотношениями Крамерса-Кронига) и будем считать, что в моменты времени г < 0 диэлектрическая проницаемость и напряжённость электрического поля равны 0. Тогда интеграл (2) можно записать следующим образом:

D (х,2,г) = |е(х,2,?)Еу (х, 2,г -?)й?.

Подставим полученное выражение (3) в уравнение (1) и воспользуемся формулой Лейбница [30]:

А Е -1

х2 у с2

де( г) . , ,,дЕу (х, 2,0) ^ '-Еу (х,2,0) + е(г)—у\ ' +

дг

дг

С Е (х, 2, г - ?) ,,

Н ? —^

г

+1.

0

дг2

= 0.

Рассмотрим планарный волновод, в который подаётся некоторый импульс. Добавим к уравнению (4) граничные и начальные условия:

'Еу (^ ^ г)| х=-4 =^1 (^ г),

Еу (^ 2 г)| 4 = ^2 (2, г) ,

2

2 е(0,I,), г £ (0, Т ]; Еу(х,2,г)|2=0 =^1(х,г),

Еу (х,2,г)12=4 =^2 (х,г), х £(- 1Х/2,1Х/2), г £(0,Т];

Еу

дг

= Ф2 (х 2),

х £

где 1х и - ширина и длина волновода соответственно; Т - время моделирования; ^ (2, г) и (2, г), ^ (х, г) и у2 (х, г), ф! (х, 2) и ф2 (х, 2) - напряжённость электрического поля на гранях волновода, на входе в волновод и выходе из него, а также в начальные моменты времени соответственно. Будем считать, что грани волновода выполнены из идеально отражающего материала (^ (2, г) = (2, г) = 0 и у2 (х, г) = 0) и в начальный момент времени поле в волноводе отсутствует (ф! (х, 2) = ф2 (х, 2) = 0), а на вход волновода поступает некоторое излучение (^ (х, г) ф 0). Из уравнения (4) видно, что первое и второе слагаемые в скобках есть произведение диэлектрической проницаемости и первого, и второго начальных условий в (7). Тогда уравнение (4) упрощается:

0

0

2

г=0

1 '

Д хъЕу- 72 N' )

д2Е (х, г,' - / ) „ —^-

д'2

= 0

Уравнение (8) вместе с граничными условиями (5), (6) и начальными условиями (7) представляет собой задачу Коши, которая будет положена в основу построения явной конечно-разностной схемы. Уравнение (8) является новой формой волнового уравнения с учётом дисперсии, в котором вторая производная по времени внесена под интеграл свёртки.

Построим на равномерной сетке явную конечно-разностную схему. Для чего в уравнении (8) произведём замену частных производных разностными соотношениями, а интеграла - интегральной суммой:

1 к „

4 Хё ] Л Ек" = Л Ек +Л Е

2 ¿^ г] ' г] х г] ъ г

7 s=0

к

ъ~г] ,

г = 1,1 -1; ] = 1, J -1; к = 1, К -1,

где Е и е] - сеточные аналоги функций Еу (х, ъ,') и

£ (х, ъ, '), взятые в узле (г,], к); ЛЕг. =

Е+1 - 2Еу + ЕЕ,-1 h2

разностный оператор Лапласа; I, J, К - количество интервалов разбиения по переменным х, ъ и '; ^ = 4 /1 hz = 1ъ / J , к = Т / К - шаги аппроксимации по переменным х, ъ и '. Выделим первое и последнее слагаемые в ряду свёртки, после чего воспользуемся начальными условиями:

КЛ'ЕТ =

4=0 к-1

= 1 ] ЛЕк5 +е к] Л $° Л'ЕЕкк =

5=1 к-1

= 1 5 Л'Е?к-5 +е 0 л Е •

С учётом (9) и (10) запишем явную конечно-разностную схему, аппроксимирующую на равномерной сетке задачу Коши (5)-(8):

„2

г- I г. \ к-1 е5

л Е = тг (Л е+ЛЕ ) - Л $5

е г] 5=1 е г

г = 1, I -1; ] = 1, J -1; к = 1, К -1;

Е] = ЕЕ] = 0, ] = 1, J -1, к = 2, К;

] 1г=0 ] 1г=/ ^

Е

1]=0

= , г = 1, I -1 , к = 2, К;

ЕЕк

= 0, г = 1, I -1, к = 2, К;

ЕЕк

= 0, г = 0, I, ] = 0, J;

ЕЕ к+1 - Е

к

= 0, г = 0, I, ] = 0, J,

где - сеточный аналог функции ^ (х, '), взятый в узле (г, ]).

2. Моделирование распространения импульса в среде из кварцевого стекла с учётом дисперсии

Проведём моделирование процесса прохождения ультракороткого импульса света внутри планарного волновода из кварцевого стекла. Для учёта зависимости диэлектрической проницаемости от частоты излучения воспользуемся моделью Селлмейера [31]: е(х, ъ, ю) = ею (х, ъ) +

+

Де„ (хг)Х2

_

<Х2-Х» (х, ъ) + 05» (х, ъ),

где £„(х, ъ) - диэлектрическая проницаемость на высоких частотах; Д£т (х, ъ) - величина резонанса; Хт (х, ъ) - резонансная длина волны; 5т (х, ъ) - коэффициент демпфирования.

Параметры модели Селлмейера для кварца представлены в табл. 1 [31]:

Таблица 1. Параметры модели Селлмейера для кварцевого стекла

т Д£т мкм X» мкм 5ОТ мкм

1 0,69616630 0,068404300 0

2 0,40794260 0,11624140 0

3 0,89747940 9,8961610 0

£„ = 1

Временные отсчёты диэлектрической проницаемости могут быть получены с помощью следующей формулы [22]:

е

где

Хт = -'1т ехр ((¿Рт - О») '), От = Пг5» / К вт = (4П2Г2 - Пг5т)1/2/^т,

Ут = 4п2Г2Дет / (4п2г2-пг5т)1/2^т

(для удобства записи зависимость от переменных х и ъ опущена).

При этом с учётом свойств дельта-функции исходное уравнение (8) можно переписать в виде:

Д Е - 4

= у г2

д'2

+

д2Еу (х,ъ,' - /) ' )--

д'2

Цх. (/ ЕУ ( X, Ъ, '

а/

= 0.

Таким образом, первое уравнение схемы (11) будет иметь вид:

ЛЕ =

' г]

-I

Х 0 +е

(Л Ек +Л ЕЕк)-

\ х г] ъ г])

=1 Х а +е

5=1

Л Е -

' г]

где Х5 - сеточный аналог х» (х, ъ), е- сеточный аналог £м(х, ъ).

Зададим начальное условие (6) в виде:

^') = cos у gup (')• gdown (')•sin(Ю0') =

0

т

5=1

0 т

2

г

5

к=0

к=0

где 8ир (г) =

Ч24, г£^0,г-);

1, г £ (ги, Т ]

г £ [0 4 - гл);

Б1П

Ф* - г]

2г„

г£

( К - га, г* ]

- функции

вхождения и затухания сигнала, 4 - время нарастания амплитуды сигнала; ^ - время спада амплитуды сигнала; 4 - время подачи сигнала; ш0 - центральная частота излучения.

Зададим параметры входного излучения для моделирования: 4=1,68 фс; /¿=1,68 фс; 4=3,36 фс. На рис. 1 представлен моделируемый импульс и его спектр.

Из рис. 1 видно, что рабочей областью спектра является промежуток длин волн от 274 до 806 нм. Построим график диэлектрической проницаемости, используя модель (12) и приведённые в табл. 1 параметры.

Из рис. 2 видно, что диэлектрическая проницаемость меняется незначительно на указанном интервале. Таким образом, разница между решениями, полученными с учётом и без учёта зависимости диэлектрической проницаемости от частоты, должна быть небольшой.

1,0

0,8

0,6

0,4

0,2

0 0,5 1,0 1,5 2,0 2,5 3,0 1,фс

Рис. 1. Интенсивность и спектр входного импульса длительностью 3,36 фс в точке (х, 2) = (0, 0) мкм

2,25 £Гам.

2,10

X, мкм

0,3 0,8

Рис. 2. График диэлектрической проницаемости для кварцевого стекла на интервале длин волн 274-806 мкм

На рис. 3 а приведена мгновенная картина интенсивности, полученная после распространения ультра-

короткого импульса внутри планарного волновода из кварцевого стекла. Для моделирования были выбраны следующие параметры: п = 1,46, = 0,532 мкм, 1х = 1 мкм, /2 = 2 мкм, Т = 3,36 фс. Для демонстрации чувствительности решения к учёту дисперсии вещества были рассчитаны 3 решения: 1) с учётом полного спектра импульса (274-806 нм), 2) с учётом укороченного спектра импульса (482-582 нм) и 3) без учёта дисперсии на основе разностного решения волнового уравнения, предложенного в [32]. Программы, реализующие конечно-разностные схемы для волнового уравнения, использовали следующие шаги аппроксимации: ^ = ^0/26 мкм, ^ = ^0/53 мкм, ^ = Х0/3192 мкм. Заметим, что для точного воспроизведения фемтосе-кундного импульса дискретность по времени должна почти в 60 раз превышать дискретность по продольной координате и более чем в 100 раз - по поперечной координате. Ошибка в расчётах составила 0,04^0.

б) 0 0,1 0,2 0,3 0,4 0,5 0,6

Рис. 3. Мгновенная осевая интенсивность импульса внутри волновода из кварца в момент времени 3,36 фс, полученная по схеме

и

1

Из рис. 3а видно, что из-за слабого изменения диэлектрической проницаемости кварца в рабочем диапазоне длин волн (на рис. 2 это изменение составляет 7%) решение волнового уравнения (14) слабо зависит от дисперсии вещества. СКО между решениями уравнения (14) с учётом дисперсии и без учёта дисперсии составило 2,87%.

Для проверки правильности решения волнового уравнения (14) было проведено сравнение с решением уравнений Максвелла FDTD-методом. Решение уравнения (14) было реализовано с помощью разработанной программы в среде MATLAB, а FDTD-метод был реализован с помощью коммерческой программы FullWAVE. Сравнению подвергались результаты, полученные с учётом частотной дисперсии материала. Для FDTD-метода были выбраны следующие параметры сетки: hx = Х0/26 мкм, hz = 1^J106 мкм, ht = Ао/5320 мкм. В качестве граничных (x = ±4/2 и z = 4) условий в программе FullWAVE здесь и далее использовались PML-условия с толщиной слоя, равной При этом для корректного сравнения полученных результатов в прграмме FullWAVE идеально отражающие условия создавались искусственно: путём помещения волновода в оболочку из материалы с высоким коэффициентом преломления (n = 1000) и толщиной в 2^о (начало PML-слоя совпадает с началом оболочки), в то время как в разработанной программе, основанной на схеме (15), идеально отражающие граничные условия выполняются на уровне самой схемы. Из рис. 3 видно, что полученные решения близки. СКО для рабочего интервала (274-806 нм) составило 6,25%, а для укороченного интервала (482-582 нм) - 3,63%. Результаты моделирования представлены на рис. 36.

3. Моделирование предвестника ультракороткого импульса в среде

В случае распространения в диспергирующей среде ультракороткого светового импульса с резким фронтом можно наблюдать появление предвестников [10, 11]. Чем круче фронт, тем большая доля энергии переносится спектральными компонентами, которые распространяются практически со скоростью света. Так, в точку zФ 0 сначала в момент времени tZ придёт ПЗ, потом в момент времени tB придёт ПБ и лишь в момент времени taR приходит оптический импульс, несущий основную энергию. Формулы для расчёта времени прибытия каждого из импульсов приведены ниже:

z

tz =-,

z

tB = - ,

z

taR = ,

где с - скорость света в вакууме; и = п / ^е (X0) - фазовая скорость; и = (дк (ю)/ дю| )-1 - групповая ско-

ч

рость; к(ю) = ю^/е (ю) / с - волновое число.

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

Для того чтобы повысить амплитуду предвестников, надо увеличить степень дисперсии материала. Рассмотрим модельный материал, в котором диэлектрическая проницаемость на спектральном рабочем интервале длин волн (274-806 мкм) будет меняться более значительно, чем диэлектрическая проницаемость для кварцевого стекла (рис. 4).

2,00

0,275 0,806

Рис. 4. График смоделированной диэлектрической проницаемости на интервале длин волн 274-806 мкм

Для моделирования также воспользуемся моделью Селлмейера (12), новые параметры которой представлены в табл. 2. Построим график диэлектрической проницаемости, используя вышеуказанную модель и приведённые в табл. 2 параметры.

Таблица 2. Параметры модели Селлмейера для модельного материала

m A£m, мкм Xm, мкм 5m, мкм

1 0,69616630 0,152956665 0

2 0,40794260 0,259923672 0

3 0,89747940 22,12848889 0

£ = 1 Ода i.

Для моделирования были выбраны следующие параметры: ги = 0 фс (резкий передний фронт импульса), 4 = 0 фс (резкий задний фронт импульса), 4 = 3,36 фс, = 0,532 мкм, 1х = 1 мкм, 12 = 8 мкм, Т = 50,4 фс, К = Х0/266 мкм, hz = Х0/354 мкм, ^ = Х0/886. На рис. 5 приведена интенсивность начального импульса с резким фронтом и его спектр.

На рис. 6 приведена зависимость интенсивности импульса от времени в точке (х, 2) = (0, 7) мкм, полученная после распространения ультракороткого импульса вида (16) с резким фронтом и длительностью 3,36 фс внутри волновода из материала, диэлектрическая проницаемость которого показана на рис. 4. Ошибка в расчётах составила 0,004^. Чтобы обнаружить предвестники с частотой почти в 10 раз большей, чем частота основного импульса, потребовалось увеличить дискретность тоже в 10 раз, по сравнению с рис. 2 и 3. Из рис. 6 видно, что основному импульсу предшествуют два предвестника. Однако интенсивность этих предвестников (в особенности ПЗ) много меньше интенсивности входного импульса (в 100 раз).

В табл. 3 и 4 приведено время гтеор и грасч появления ПЗ и ПБ и основного импульса в точках (х, 2) = (0, 3) и (х, 2) = (0, 7) мкм, полученное с помощью вышеприведённых формул (17), а также рассчитанное на основе данных моделирования.

\Е\2, а. и.

0 12 3

Рис. 5. Интенсивность и спектр входного импульса с резким фронтом длительностью 3,36 фс в точке (х, 1) = (0, 0) мкм

8,5 9,0 9,5 10,0 10,5 11,0 11,5

Рис. 6. Временная зависимость интенсивности импульса в точке (х, 1) = (0, 7) мкм

Из табл. 3 и 4 видно, что отклонение для ПБ достаточно мало, а для ПЗ достаточно велико. При этом чем дальше от входа в волновод находится область наблюдения, тем ниже ошибка полученных результатов для предвестников (отклонение уменьшается практически в 2 раза). Различие между отклонением, полученным в случае расчёта времени появления предвестника по началу импульса, и отклонением, полученным в случае расчёта времени появления предвестника по центру импульса, составляет менее 2%. Это означает, что ошибка (отклонение) не зависит от точности определения начала предвестника.

Из рис. 5 видно, что задний фронт импульса обрывается не на нулевом значении интенсивности. Возникло предположение, что именно резкий задний фронт приводит к высокой задержке по времени предвестника Зоммерфельда и к большому отклонению рассчитанного времени задержки от теоретического. В связи с этим было проведено дополнительное моделирование, при котором длина импульса была увеличена с 3,34 фс до 3,55 фс, чтобы в импульсе уложилось ровно две длины волны, а задний фронт обрывался на нулевом значении интенсивности. В результате модели-

рования также наблюдался предвестник Зоммерфель-да. При этом отклонение времени прибытия предвестника, полученного на основе моделирования, от времени, рассчитанного с помощью формул (17), для точки (х, 1) = (0, 3) мкм снизилось почти в 8 раз и составило 5,1%, а для точки (х, 1) = (0, 7) мкм снизилось почти в 4 раза и составило 5,14%.

Таблица 3. Время появления предвестников и основного импульса в точке (х, 1) = (0, 3) мкм

Тип импульса 1теор, фс 1эксп, фс Отклонение, %

Предвестник Зоммерфельда 10,00 13,9 39

Предвестник Бриллюэна 15,16 14,375 5,19

Несущий импульс 16,72 16,64 0,49

Таблица 4. Время появления предвестников и основного импульса в точке (х, 1) = (0, 7) мкм

Тип импульса 1теор, фс 1эксп, фс Отклонение, %

Предвестник Зоммерфельда 23,35 28,02 20

Предвестник Бриллюэна 35,37 34,326 2,95

Несущий импульс 39,02 37,06 5,03

Заключение

В работе получены следующие результаты:

- волновое уравнение размерностью (2+1) записано в нестандартной форме: вторая производная по времени от напряжённости электрического поля находится под знаком интеграла свёртки, описывающим временную зависимость электрической индукции (уравнение (8));

- предложена явная конечно-разностная схема (11) для решения волнового уравнения (8), учитывающего зависимость диэлектрической проницаемости от частоты; алгоритм решения уравнения (8) реализован в виде программы в среде МА^АВ;

- проведено сравнение результатов моделирования процесса распространения фемтосекундного импульса длительностью 3,36 фс с несущей длиной волны 532 нм (ширина спектра 274-806 нм) в планарном волноводе из кварцевого стекла, дисперсия которого описывается моделью Селлмейера, полученных с учётом и без учёта дисперсии; СКО составило 3% (рис. 3а);

- проведено также сравнение результатов численного моделирования, распространения того же импульса в той же среде, полученных с помощью решения волнового уравнения (8) и решения уравнений Максвелла методом FDTD с учётом дисперсии, реализованном в программе FullWAVE; СКО составило 6% (рис. 3б);

- моделирование распространения фемтосекунд-ного импульса с резким передним фронтом в планар-ном волноводе из материала, частотная дисперсия которого описывается на базе модели Селлмейера, по-

казало наличие разделённых во времени предвестников Зоммерфельда и Бриллюэна, интенсивность которых в 100 раз меньше интенсивности основного импульса; СКО времени появления предвестников, полученного с помощью теоретических формул, от времени, полученного при моделировании, составило для ПЗ и ПБ 20% и 3 %.

Благодарности

Работа выполнена при поддержке Министерства образования и науки Российской Федерации, гранта Президента РФ поддержки ведущих научных школ (НШ-4128.2012.2) и гранта РФФИ (13-07-97008).

Авторы выражают благодарность рецензенту за высказанную гипотезу о влиянии заднего фронта импульса на время прибытия предвестника.

Литература

1. Hecht, J. Spectral Broadening Advances Quest for Single-Cycle Pulses / J. Hecht // Laser Focus World. - 2011. -V. 47(8). - P. 65-70.

2. Fourmaux, S. Laser Pulse Contrast Ratio Cleaning in 100 TW Scale Ti: Sapphire Laser Systems / S. Fourmaux, S. Payeur, Ph. Lassonde, J.C. Kieffer and F. Martin // Laser Systems for Applications. - 2011. - P. 139-154.

3. Block, М. Few-cycle high-contrast vortex pulses / M. Block, J. Jahns and R. Grunwald // Optics Letters. -2012. - V. 37(18). - P. 3804-3806.

4. Povolotskiy, A. 2D and 3D laser writing for integrated optical elements creation / A. Povolotskiy, A. Shimko, A. Man-shina // Mondello Proceedings of WF0PC2005: Elioticinese Service Point srl. - 2005. - V. 4. - P. 196-202.

5. Cheng, Ya. Microfabrication of 3D hollow structures embedded in glass by femtosecond laser for Lab-on-a-chip applications / Ya. Cheng, K. Sugioka, K. Midorikawa // Applied Surface Science. - 2005. - V. 248 - P. 172-176.

6. Liu, X. Laser ablation and micromachining with ultrashort laser pulses / X. Liu, D. Du, G. Mourou // IEEE Quantum Electron. - 1997. - V. 38. - P. 1706.

7. Крюков, П.Г. Лазеры ультракоротких импульсов и их применения / П.Г. Крюков. - Долгопрудный: Издательский дом «Интеллект», 2012. - 248 с.

8. Бочкарёв, Н.Н. Взаимодействие фемтосекундных лазерных импульсов с биологическим веществом / Н.Н. Бочкарёв [и др.]. - Томск: ТПУ, 2007. - 121 с.

9. Liu, Z. Ultraviolet conical emission produced by highpower femtosecond laser pulse in transparent media / Z. Liu, X. Lu, Q. Liu, S. Sun, L. Li, X. Liu, B. Ding, B. Hu // Appl. Phys. B. - 2012. - V. 108 - P. 493-500.

10. Piglosiewicz, B. Ultrasmall bullets of light-focusing few-cycle light pulses to the diffraction limit / B. Piglosiewicz, D. Sadiq, M. Mascheck, S. Schmidt, M. Silies, P. Vasa and C. Lienau // Optics Express. - 2011. - V. 19(15) -P. 14451-14463.

11. Ахманов, C.A. Оптика фемтосекундных лазерных импульсов / С.А. Ахманов, В.А. Выслоух, А.С. Чиркин. -М.: Наука, 1988. - 312 с.

12. Oughstun, K.E. Electromagnetic pulse propagation in cau-sial dielectrics / K.E. Oughstun, G.C. Sherman. - SpringerVerlag, 1994. - 465 p.

13. Li, C. Approach to accurately measuring the speed of optical precursors / C. Li, Z. Zhou, H. Jeong, G. Guo // Phys. Rev. A. - 2011. - V. 84 - P. 043803.

14. Safian, R. Joint time-frequency and FDTD analysis of precursor fields in dispersive media / R. Safian, C.D. Sarris, M. Mojahedi // Phys. Rev. E. - 2006.- V. 73.- P. 066602.

15. Jeong, H. Evolution of Sommerfeld and Brillouin precursors in intermediate spectral regimes / H. Jeong, U.L. Osterberg, T. Hansson // J. Opt. Soc. Am. B. - 2009. - V. 26. -P. 2455-2460.

16. Macke, B. From Sommerfeld and Brillouin forerunners to optical precursors / B. Macke and B. Segard // Phys. Rev. A. - 2013. - V. 87 - P. 043830.

17. Macke, B. Simle asymptotic forms for Sommerfelod and Brillouin precursore / B. Macke and B. Segard // Phys. Rev.

A. - 2012. - V. 86 - P. 013837.

18. Luebbers, R.J. A frequency-depended finite-difference time-domain formulation for dispersive materials / R.J. Luebbers, F.P. Hunsberger, K.S. Kunz, R.B. Standler and M. Schneider // IEEE Trans. Electromagn. Compat. -

1990. -V. 32. - P. 222-227.

19. Luebbers, R.J. FDTD for Nth-order dispersive media / R.J. Luebbers and F.P. Hunsberger // IEEE Trans. Antennas. Propagat. - 1992. -V. 40. - P. 1297-1301.

20. Bui, M.D. Propagation of transients in dispersive dielectric media / M.D. Bui, S.S. Stuchly and G.I. Costache // IEEE Trans. Microwve Theory Technol. - 1991. - V. 39 -P. 1165-1171.

21. Hawkins, R.J. Linear electronic dispersion and finite-difference time-domain calculations: a simple approach / R.J. Hawkins and J.S. Kallman // J. Lightwave Technol. -1993. - V. 11 - P. 1872-1874.

22. Kelley, D.F. Piecewise linear recursive convolution for dispersive media using FDTD / D.F. Kelley and R.J. Luebbers -IEEE Trans. Antennas. Propagat. - 1996. -V. 44. - P. 792-797.

23. Kashiwa, T. A treatment by the FD-TD method of the dispersive characteristics associated with electronic polarizaion / T. Kashiwa and I. Fukai // Microwave Opt. Technol. Lett. - 1990. - V. 3(5) - P. 203-205.

24. Kashiwa, T. A finite-difference time-domain formulation for transient propagation in dispersive media associated with Cole-Cole's circular arc law / T. Kashiwa, Y. Ohtomo and I. Fukai // Microwave Opt. Technol. Lett. - 1990. -V. 3(12) - P. 16-419.

25. Joseph, R.M. Direct time integration of Maxwell's equations in linear dispersive media with absorption for scattering and propagation of femtosecond electromagnetic pulses / R.M. Joseph, S.C. Hagness and A. Taflove // Opt. Lett. -

1991. - V. 16. - P. 1412-1414.

26. Gandhi, O.P. A frequency-depended finite-difference timedomain formulation for general dispersive media / O.P. Gandhi, B.-Q. Gao and J.Y. Chen // IEEE Trans. Mi-crowve Theory Technol. - 1993. - V. 41. - P. 658-664.

27. Korner, T.O. Auxiliary differential equation: efficient implementation in the finite-difference time-domain method / T.O. Korner and W. Fichtner // Opt. Lett. - 1997. -V. 22(21) - P. 1586-1588.

28. Liu, Y. Formulation of the finite-difference time-domain method for the analysis of axially symmetric metal nanode-vices / Y. Liu and W. Yu // Journal of Modern Optics. -2012. - V. 59(16) - P. 1439-1447.

29. Sullivan, D.M. Frequency-depended FDTD methods using Z transform / D.M. Sullivan // IEEE Trans. Antennas. Propagat. - 1992. -V. 40. - P. 1223-1230.

30. Ильин, В.А. Математический анализ. Ч.1. / В.А. Ильин,

B.А. Садовничий, Бл.Х. Сендов. - М: Изд-во МГУ, 1985. - 663 с.

31. Couairon, A. Filamentation and damage in fused silica induced by tightly focused femtosecond laser pulses / A. Cou-

airon, L. Sudrie, M. Franco, B. Prade, A. Mysyrowicz // Phys. Rev. B. - 2005. - V. 71. - P. 125435-125441.

32. Козлова, Е.С. Моделирование распространения короткого двумерного импульса света / Е.С. Козлова, В.В. Котляр // Компьютерная оптика. - 2012. - Т. 36, № 2. - С. 158-164.

References

1. Hecht, J. Spectral Broadening Advances Quest for Single-Cycle Pulses / J. Hecht // Laser Focus World. - 2011. -V. 47(8). - P. 65-70.

2. Fourmaux, S. Laser Pulse Contrast Ratio Cleaning in 100 TW Scale Ti: Sapphire Laser Systems / S. Fourmaux, S. Payeur, Ph. Lassonde, J.C. Kieffer and F. Martin // Laser Systems for Applications. - 2011. - P. 139-154.

3. Block, М. Few-cycle high-contrast vortex pulses / M. Block, J. Jahns and R. Grunwald // Optics Letters. -2012. - V. 37(18). - P. 3804-3806.

4. Povolotskiy, A. 2D and 3D laser writing for integrated optical elements creation / A. Povolotskiy, A. Shimko, A. Man-shina // Mondello Proceedings of WF0PC2005: Elioticinese Service Point srl. - 2005. - V. 4. - P. 196-202.

5. Cheng, Ya. Microfabrication of 3D hollow structures embedded in glass by femtosecond laser for Lab-on-a-chip applications / Ya. Cheng, K. Sugioka, K. Midorikawa // Applied Surface Science. - 2005. - V. 248 - P. 172-176.

6. Liu, X. Laser ablation and micromachining with ultrashort laser pulses / X. Liu, D. Du, G. Mourou // IEEE Quantum Electron. - 1997. - V. 38. - P. 1706.

7. Krukov, P.G. Lasers of ultrafast pulses and it's application / P.G. Krukov. - Dolgoprudniy: "Intellekt" Publishing House, 2012. - 248 p. - (In Russian).

8. Bohkarev, N.N. The interaction of femtosecond laser pulses with biological substance / N.N. Bockarev [et al.]. -Tomsk: "TPU" Publishing, 2007. - 121 p. - (In Russian).

9. Liu, Z. Ultraviolet conical emission produced by highpower femtosecond laser pulse in transparent media / Z. Liu, X. Lu, Q. Liu, S. Sun, L. Li, X. Liu, B. Ding, B. Hu // Appl. Phys. B. - 2012. - V. 108 - P. 493-500.

10. Piglosiewicz, B. Ultrasmall bullets of light-focusing few-cycle light pulses to the diffraction limit / B. Piglosiewicz, D. Sa-diq, M. Mascheck, S. Schmidt, M. Silies, P. Vasa and C. Lie-nau // Optics Express. - 2011. - V. 19(15) - P. 14451-14463.

11. Ahmanov, S.A. Optics of Femtosecond Laser Pulses / S.A. Ahmanov, V.A. Vislouh, A.S. Chirkin. - Moscow: "Nauka" Publisher, 1988. - 312 p. - (In Russian).

12. Oughstun, K.E. Electromagnetic pulse propagation in cau-sial dielectrics / K.E. Oughstun, G.C. Sherman. - SpringerVerlag, 1994. - 465 p.

13. Li, C. Approach to accurately measuring the speed of optical precursors / C. Li, Z. Zhou, H. Jeong, G. Guo // Phys. Rev. A. - 2011. - V. 84 - P. 043803.

14. Safian, R. Joint time-frequency and FDTD analysis of precursor fields in dispersive media / R. Safian, C.D. Sarris, M. Mojahedi // Phys. Rev. E. - 2006.- V. 73.- P. 066602.

15. Jeong, H. Evolution of Sommerfeld and Brillouin precursors in intermediate spectral regimes / H. Jeong, U.L. Osterberg, T. Hansson // J. Opt. Soc. Am. B. - 2009. - V. 26. -P. 2455-2460.

16. Macke, B. From Sommerfeld and Brillouin forerunners to optical precursors / B. Macke and B. Segard // Phys. Rev. A. - 2013. - V. 87 - P. 043830.

17. Macke, B. Simle asymptotic forms for Sommerfelod and Brillouin precursore / B. Macke and B. Segard // Phys. Rev. A. - 2012. - V. 86 - P. 013837.

18. Luebbers, R.J. A frequency-depended finite-difference time-domain formulation for dispersive materials / R.J. Luebbers, F.P. Hunsberger, K.S. Kunz, R.B. Standler and M. Schneider // IEEE Trans. Electromagn. Compat. -

1990. -V. 32. - P. 222-227.

19. Luebbers, R.J. FDTD for Nth-order dispersive media / R.J. Luebbers and F.P. Hunsberger // IEEE Trans. Antennas. Propagat. - 1992. -V. 40. - P. 1297-1301.

20. Bui, M.D. Propagation of transients in dispersive dielectric media / M.D. Bui, S.S. Stuchly and G.I. Costache // IEEE Trans. Microwve Theory Technol. - 1991. - V. 39 - P. 1165-1171.

21. Hawkins, R.J. Linear electronic dispersion and finite-difference time-domain calculations: a simple approach / R.J. Hawkins and J.S. Kallman // J. Lightwave Technol. -1993. - V. 11 - P. 1872-1874.

22. Kelley, D.F. Piecewise linear recursive convolution for dispersive media using FDTD / D.F. Kelley and R.J. Luebbers -IEEE Trans. Antennas. Propagat. - 1996. -V. 44. - P. 792-797.

23. Kashiwa, T. A treatment by the FD-TD method of the dispersive characteristics associated with electronic polarizaion / T. Kashiwa and I. Fukai // Microwave Opt. Technol. Lett. - 1990. - V. 3(5) - P. 203-205.

24. Kashiwa, T. A finite-difference time-domain formulation for transient propagation in dispersive media associated with Cole-Cole's circular arc law / T. Kashiwa, Y. Ohtomo and I. Fukai // Microwave Opt. Technol. Lett. - 1990. -V. 3(12) - P. 16-419.

25. Joseph, R.M. Direct time integration of Maxwell's equations in linear dispersive media with absorption for scattering and propagation of femtosecond electromagnetic pulses / R.M. Joseph, S.C. Hagness and A. Taflove // Opt. Lett. -

1991. - V. 16. - P. 1412-1414.

26. Gandhi, O.P. A frequency-depended finite-difference timedomain formulation for general dispersive media / O.P. Gandhi, B.-Q. Gao and J.Y. Chen // IEEE Trans. Microwve Theory Technol. - 1993. - V. 41. - P. 658-664.

27. Korner, T.O. Auxiliary differential equation: efficient implementation in the finite-difference time-domain method / T.O. Korner and W. Fichtner // Opt. Lett. - 1997. -V. 22(21) - P. 1586-1588.

28. Liu, Y. Formulation of the finite-difference time-domain method for the analysis of axially symmetric metal nanode-vices / Y. Liu and W. Yu // Journal of Modern Optics. -2012. - V. 59(16) - P. 1439-1447.

29. Sullivan, D.M. Frequency-depended FDTD methods using Z transform / D.M. Sullivan // IEEE Trans. Antennas. Propagat. - 1992. -V. 40. - P. 1223-1230.

30. Ilin, V.A. Mathematical analysis. P.1. / V.A. Ilin, V.A. Sa-dovnishii, Bl.H. Sendov. - Moscow: "MSU Press" Publisher, 1985. - 663 p. - (In Russian).

31. Couairon, A. Filamentation and damage in fused silica induced by tightly focused femtosecond laser pulses / A. Cou-airon, L. Sudrie, M. Franco, B. Prade, A. Mysyrowicz // Phys. Rev. B. - 2005. - V. 71. - P. 125435-125441.

32. Kozlova, E.S. Simulation of ultrafast 2D light pulse / E.S. Kozlova, V.V. Kotlyar // Computer Optics. - 2012. -V. 36(2) - P. 158-164. - (In Russian).

SIMMULATIONS OF SOMMERFELD AND BRILLOUIN PRECURSORS IN THE MEDIUM WITH FREQUENCY DISPERSION USING NUMERICAL METHOD OF SOLVING WAVE EQUATIONS

E.S. Kozlova, V.V. Kotlyar Image Processing Systems Institute of the RAS

Abstract

The wave equation (2 +1) taking into account the frequency dispersion of the material for TE polarization is written in a form in which the second time derivative included under the sign of the convolution integral describing the electric induction. An algorithm for solving such equations using explicit finite difference scheme is developed. The solution of the wave equation, describing the propagation of ultrashort pulse duration of 3.36 fs (in the spectral range 274-806 nm) in a planar waveguide made of quartz glass at the center wavelength of 532 nm, is different from solutions of this equation without dispersion for 3% and from solutions of Maxwell's equations derived by FDTD-method in the program FullWAVE, for 6%. A detailed analysis of the obtained solutions for a femtosecond pulse with a sharp initial front resulted in the discovery of precursors, which come to the observation point before the main pulse and the intensity is 100 times smaller. Moreover, the calculated time delay Sommerfeld precursor differs from the theoretical to 20%, And the Brillouin precursor of 3%.

Key words: wave equation, permettivety, explicit finite-difference scheme, numerical simulation, ultrashort pulse, the optical precursor.

Сведения об авторах

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

E-mail: kozlova. elena. s@gmail. com .

Elena Sergeevna Kozlova, master of mathematics and computer science. Currently studies at Samara State Aerospace University. Research interests are diffractive optics and numerical methods.

Котляр Виктор Викторович, 1957 года рождения, доктор физико-математических наук, профессор, заведующий лабораторией лазерных измерений ИСОИ РАН и, по совместительству, профессор кафедры технической кибернетики Самарского государственного аэрокосмического университета (СГАУ). В 1979 году окончил физический факультет Куйбышевского государственного университета, в 1988 году защитил кандидатскую диссертацию в Саратовском государственном университете, а в 1992 году - докторскую диссертацию в Центральном конструкторском бюро Уникального приборостроения РАН (г. Москва). Область научных интересов: нанофотоника, дифракционная компьютерная оптика. Публикации: 300 научных трудов, 5 монографий, 7 авторских свидетельств.

Victor Victorovich Kotlyar is a head of Laboratory at the Image Processing Systems Institute (Samara) of the Russian Academy of Sciences and professor of Computer Science Department at Samara State Aerospace University. He received his MS, PhD and DrSc degrees in physics and mathematics from Samara State University (1979), Saratov State University (1988) and Moscow Central Design Institute of Unique Instrumentation, the Russian Academy of Sciences (1992). He is SPIE and OSA member. He is co-author of 300 scientific papers, 5 books and 7 inventions. His current interests are diffractive optics, gradient optics, nanophotonics, optical vortices.

Поступила в редакцию 8 марта 2013 г.