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

УДК 535.3

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

© 2013 Е.С. Козлова1- 2

1 Институт систем обработки изображений РАН, г.Самара 2 Самарский государственный аэрокосмический университет имени академика С.П. Королева (национальный исследовательский университет)

Поступила в редакцию 23.11.2013

Предложен алгоритм решения волнового уравнения для ТЕ-поляризации с учётом частотной дисперсии материала, с помощью явной разностной схемы. Анализ численного моделирования распространения прямоугольного фемтосекундного импульса в кварце позволил обнаружить предвестники, которые приходят в точку наблюдения раньше основного импульса в 1,4 раза и по интенсивности в 700 раз меньше его. Рассчитанное время задержки предвестников Зоммерфельда и Бриллюэна отличается от теоретического на 7%.

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

1. ВВЕДЕНИЕ

Ультракороткие импульсы имеют широкое применение в современной науке и технике [1 ]. Для моделирования распространения света в дисперсных материалах были разработаны различные модификации БЭТО-метода, позволяющие в процессе вычислений учитывать зависимость диэлектрической проницаемости от частоты [2].

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

Т7к ц! +1/2

екций поля: Еу , Нх

Hk +1/2 . Проведено

фронтом в планарном волноводе из кварцевого стекла с учётом дисперсии показало наличие предвестников Зоммерфельда и Бриддюэна, разделённых в пространстве.

2. КОНЕЧНО-РАЗНОСТНАЯ СХЕМА ДЛЯ РЕШЕНИЯ ВОЛНОВОГО УРАВНЕНИЯ С УЧЕТОМ ДИСПЕРСИИ

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

1 э 2Д, „

(1)

A E -

xz y

Э t2

= 0

э2 Э2

где A xz = +

Эх2 Эz2

оператор Лапласа, E

сравнение результатов численного моделирования, полученных с помощью разработанной схемы и программного пакета FullWAVE (RSoft), реализующего дисперсионный FDTD-метод. Моделирование распространения фемтосекундного прямоугольного импульса с резким передним

Козлова Елена Срегеевна, магистр прикладной математики и информатики, аспирант кафедры технической кибернетики, стажер-исследователь. E-mail: kozlova.elena.s@gmail.com

проекция вектора напряжённости электрического поля на ось у (ТЕ-поляризация); x и 2 - поперечная и продольная пространственные координаты; Ь - время; с - скорость света в вакууме; Бу - проекция вектора электрической индукции на ось у. В случае немагнитной среды электрическая индукция будет представлена интегралом свёртки [3]:

Ву (х, 2, г) = \е (х, 2, г) Еу (х, 2, г - г) £. (2)

0

Будем считать, что свет распространяется в па-ланарном волноводе, грани которого выполнены из идеально отражающего материала, и в начальный момент времени поле в волноводе отсутствует.

Воспользуемся принципом причинности и будем считать, что в моменты времени Ь<0 диэ-

лектрическая проницаемость и напряженность электрического поля равны 0:

г

Бу (х, 2,г) = | '(х, 2, г )Еу (х, г,г - г )£. (3)

0

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

1 \ ,д2Е (х,;

х, 2, г - г

= 0. (4)

Для учета зависимости диэлектрической проницаемости от частоты излучения воспользуемся следующей обобщенной моделью:

' (^ г, ю)=£^(x, г )+ Ё хт (х, ^ ®), (5)

т

где £ (х, г) - диэлектрическая проницаемость на высоких частотах; Хт (х, г, а>) - резонансная составляющая. Данная модель, путем определения формы зависимости резонансной составляющей от частоты, преобразуется во все общепринятые модели дисперсии: Лоренца, Дебая, Дру-де, Селмейера.

Временные отсчеты диэлектрической проницаемости могут быть получены с помощью обратного преобразования Фурье:

' (х, г,г) = ' (х, г)8 (г)+ Ё %т (х, 2,г), (6)

т

где 8 (г) - дельта-функция.

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

—| £ л ,Ек + Ё х- л ,Ек ^1 | = л Ек + л Ек ,(7)

2 I Ч г Ч /—I лЧ г Ч I х Ч г Ч 'V '

С V * = 0 У

где Ек, £' , - сеточные аналоги функций

Б(х,г,Ь), £ (х, г, г) и хт (х, г, г), взятые в узле

г Е.,, - 2Е . + Е.

(г,3,к); ЛЕ. -

'/г'+1 г г-1

И2

разностный

В итоге мы получили явную конечно разностную схему для волнового уравнения с учетом зависимости диэлектрической проницаемости от частоты:

ЛЕ =

г .

■(Л Е+л Е)-Ё *. (9)

Х- + '■■ ] ]

=, Х- +'-

2. МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ИМПУЛЬСА В ВОЛНОВОДЕ ИЗ КВАРЦЕВОГО СТЕКЛА

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

А£ (х г)л

(х, г)

-,(10)

£(х,г,Л)=£ (х,г)+ Ё--, о, ч

1" ' тл-л(х,г)+аЛтх. ,

где Л - длинна волны; £„ (х, г) - диэлектрическая проницаемость на высоких частотах; А'т (х, г) - величина резонанса; Лт (х, г) - резонансная длина волны; I - мнимая единица; Г)т (х, г) - коэффициент демпфирования. Параметры модели Селлмейера для кварца взяты из работы [4].

Резонансная составляющая в (6) с учетом (10) представляет собой экспоненциальную функцию (для удобства записи зависимость от переменных х и 2 опущена) (Цт (х, г) =0):

Хт (г) =-1

2лс А£

_т

Л

-ехр

( 2тсг ^

(11)

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

¥(х, г) =

( тх ^

008

• 81П

( т ^

V у

■8Ш (в)ог), (12)

2тс

оператор Лапласа; к - шаг аппроксимации; I, ], К - количество интервалов разбиения по переменным х, 2 и V; И = —, И = — и Иг = —г--шаги

х I г 3 г К

аппроксимации по переменным х, 2 и V; 1х и ¡г -ширина и длина волновода; Т - время моделирования. Выделим первое и последнее слагаемые в ряду свертки, после чего воспользуемся начальными условиями:

Ё Х]лЕк * = Ё Х]лЛк] * + %%]лЕ . (8)

где V - время подачи импульса; (О0 =--цент' Л0 ральная частота излучения; Л0 - центральная длина волны излучения.

Зададим параметры входного излучения для моделирования: Ь=3,55 фс. На рис. 1а представлен начальный импульс и его спектр. Для проверки правильности решения волнового уравнения (4) было проведено сравнение с решением уравнений Максвелла РОТЭ-методом. Решение уравнения (4) было реализовано с помощью разработанной программы в среде МЛТЬЛВ, а РОТЭ-метод был реализован с помощью коммерческой программы FullWЛVE. Сравнению подвергались результаты, полученные с учетом частотной дисперсии материала. В программе, реализующая

конечно-разностную схему для волнового уравнения, использовались следующие шаги аппроксимации: йх= Л0 /26 мкм, Ь,= Л0 /53 мкм, й£=Л0 /3192 мкм. Для БЭТЭ-метода были выбраны следующие параметры сетки: Ь,= Л0/26 мкм, Н= Л0 /106 мкм, й=Л0 /5320 мкм. Из рисунка 1б видно, что полученные решения близки. Среднеквадратичес-кое отклонение составило 3,6%.

3. МОДЕЛИРОВАНИЕ ПРЕДВЕСТНИКОВ ИМПУЛЬСА С РЕЗКИМ ПЕРЕДНИМ ФРОНТОМ

В случае распространения в диспергирующей среде ультракороткого светового импульса с резким передним фронтом можно наблюдать появление предвестников [3]. Так в точку 2 ^ 0 сначала в момент времени Ьг придёт предвестник Зоммерфельда, потом в момент Ьв - предвестник Бриллюэна и лишь в момент Ьск - оптический импульс. Формулы для расчёта времени прибытия каждого из импульсов даны в [3]:

_ 2 _ 2 _ 2

= ; гв = ; гок = , (12)

С V и

где V - фазовая скорость; и - групповая скорость.

Промоделируем распространение прямоугольного импульса в волноводе из кварцевого стекла. Зададим начальное условие в виде:

\Е\2,а.и.

¥ (х, г) =

008

( ТГ^Г \

тех

т~

• гесг

Г г >

V /

• 81П

К г ).(13)

Для моделирования были выбраны следующие параметры:, Ь=3,55 фс, Л0 =0,532 мкм, 1=1 мкм, 1=30 мкм, Т=100,1 фс, Н= Л0/266 мкм, Н= Л0 /354 мкм, Н= Л0/886. На рис. 2а приведена интенсивность начального импульса с резким передним и задним фронтом и его спектр. На рисунке 2б приведена зависимость интенсивности импульса от времени в точке (х,2) = (0,7) мкм, полученная после распространения ультракороткого импульса вида (13) длительностью 3,55 фс внутри волновода из кварцевого стекла.

Ошибка в расчётах составила - 0,004 Л0.

Пройдя всего 7 мкм в волноводе, начальный импульс (рис. 2а) и его начальный спектр кардинально изменились. Из рисунок 2б видно, что основному импульсу предшествуют два предвестника. Однако интенсивность этих предвестников (в особенности предвестник Зоммерфельда) много меньше интенсивности несущего импульса (в 700 раз). На рис. 2в приведен спектр импульса, показанного на рисунке 2б. Из рис. 2в видно, что в спектре преобладают низкие частоты, характерные для предвестника Бриллюэна, в то время как высокие частоты, характерные для предвестника Зоммерфельда, гораздо слабее.

\Е\,а.и.

Рис. 1. а - интенсивность |Еу| и амплитуда спектра 8 входного импульса в точке (х^)=(0,0) мкм; б - мгновенная интенсивность импульса вдоль оси z, полученная: по схеме (9) (линия 1) и РСТС-методом (линия 2)

1,9

\Е\, а.и.

£ а.и.

3,55

13\——1Г

.¿м

7 а, мкм ю,5 ^_

мкм

б)

Рис. 2. а - интенсивность Е и амплитуда спектра 8 входного прямоугольного импульса в точке (х^)=(0,0) мкм;

б - временная зависимость интенсивности Е ; в - спектр 8 импульса в точке (х^) = (0,7) мкм

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

Тип импульса ^теор, фс ^расч-> фс Отклонение, %

Предвестник Зоммерфельда 23,35 25,09 7,45

Предвестник Бриллюэна 35,37 32,89 7,53

Несущий импульс 39,02 34,19 12,38

Несущая частота С00 так же не вносит основной вклад в спектр импульса. В табл. 1 приведено время V и V появления предвестников Зом-

Г теор расч г

мерфельда и Бриллюэна и основного импульса в точке (х,2) = (0,7) мкм, полученные с помощью формул (12), а также рассчитанные на основе данных моделирования.

Из табл. 1 видно, что предвестник Зоммер-фельда приходит точку наблюдения в 1,36 раз быстрее, чем основной импульс. эта цифра согласуется с показателем преломления кварца на несущей длине волны п0=1,46 ( Л0 =532 нм).

4. ЗАКЛЮЧЕНИЕ

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

2) предложена явная конечно-разностная схема (9) для решения волнового уравнения (4), учитывающего зависимость диэлектрической проницаемости от частоты; алгоритм решения уравнения (4) реализован в виде программы в среде МЛТЬЛВ;

3) проведено сравнение результатов численного моделирования, распространения фемтосекун-дного импульса длительностью 3,55 фс с длиной волны 532 нм в планарном волноводе из кварцевого стекла, дисперсия которого описывается моделью Селлмейера, полученных с помощью решения

волнового уравнения (4), и решения уравнений Максвелла методом FDTD с учётом дисперсии, реализованном в программе FullWAVE; среднеквад-ратическое отклонение составило 4% (рис. 1б);

4) моделирование распространения фемтосе-кундного прямоугольного импульса длительностью 3,55 фс в том же волноводе, показало наличие разделенных во времени предвестников Зом-мерфельда и Бриллюэна, интенсивность которых в 700 раз меньше интенсивности несущего импульса (рис. 2б); среднеквадратическое отклонение времени появления предвестников, полученного с помощью теоретических формул, от времени, полученного при моделировании, составило для предвестника Зоммерфельда и предвестника Бриллюэна около 7% (табл. 1).

СПИСОК ЛИТЕРАТУРЫ

1. Few-cycle high-contrast vortex pulses / M. Block, J. Jahns, and R. Grunwald // Optics Letters. 2012. Vol. 37, № 18. P. 3804-3806.

2. Liu Y., Yu W. Formulation of the finite-difference timedomain method for the analysis of axially symmetric metal nanodevices // Journal of Modern Optics. 2012. Vol. 59, №16. P. 1439-1447.

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

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

SIMULATION OF FEMTOSECOND PULSES IN THE MEDIUM WITH FREQUENCY DISPERSION

© 2013 E.S. Koziova1,2

1 Image Processing Systems Institute RAS, Samara 2 Samara State Aerospace University named after Academician S.P. Koroiyov (National Research University)

An algorithm for solving the wave equation for the TE polarization, which is taking into account the frequency dispersion of the material and using the explicit difference scheme, is proposed. Analysis of the numerical simulation of propagation of femtosecond rectangular pulse in a quartz allowed to find precursors that come to the observation point for 1.4 times faster than the main pulse and the intensity of this precursors is in 700 times smaller than the intensity of the main pulse. The calculated delay time of Sommerfeld and Brillouin precursors differents from the theoretical to 7%.

Key words: mathematical model, wave equation, frequency dispersion, Sellmeyer's model, femtosecond rectangular pulse.

Elena Kozlova, Master of Mathematics and Computer Science, Graduate Student, Trainee Researcher. E-mail: kozlova.elena.s@gmail.com