CC BY
18
ТЕОРЕТИЧЕСКАЯ И МАТЕМАТИЧЕСКАЯ ФИЗИКА
Исследование распространения электромагнитных импульсов через фотонные кристаллические структуры
А. Н. Боголюбов0, И. А. Буткарев, Ю. С. Дементьева
Московский государственный университет имени М. В. Ломоносова, физический факультет, кафедра математики. Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2. E-mail: а bogan7@yandex.ru
Статья поступила 08.07.2010, подписана в печать 30.08.2010
Построена и применена математическая модель исследования двумерных и трехмерных фотонных кристаллов и основанных на них волноведущих систем. Численный алгоритм построения решения основан на таких численных методах, как метод конечных разностей во временной области (FDTD), метод TF/SF и метод идеально согласованного слоя (PML), что позволяет применять данный алгоритм к исследованию и других волноведущих систем. Приводятся результаты моделирования конкретных волноведущих систем.
Ключевые слова: фотонные кристаллы, математическое моделирование, параллельные вычисления. УДК: 681.7.068. PACS: 02.60.Cb.
Ведение
В настоящей работе строится и применяется математическая модель исследования волноведущих систем на основе фотонных кристаллов. Концепция фотонных кристаллов (ФК) была впервые предложена в 1987 г. Э. Яблоновичем [1, 2]. В общем случае фотонный кристалл — это материал, структура которого характеризуется периодическим изменением коэффициента преломления. Как известно, кристаллы всех типов могут эффективно рассеивать некоторые виды излучения при условии, что параметры решетки кристалла имеют тот же порядок, что и длина волны излучения. Аналогичным образом, будучи прозрачными для широкого спектра электромагнитного излучения, фотонные кристаллы не пропускают свет с длиной волны, сравнимой с периодом структуры фотонного кристалла. Эти спектральные диапазоны получили название «фотонные запрещенные зоны» (photonic bandgap, PBG). Если запрещенные зоны для всех без исключения направлений в кристалле перекрываются, то можно говорить о полной запрещенной зоне. Наличие полной запрещенной зоны означает, что в определенном спектральном диапазоне свет любой поляризации не может распространяться в ФК ни в каком направлении в одном, двух или трех измерениях.
Численный алгоритм
Распространение электромагнитного поля в фотон-но-кристаллической системе описывается уравнениями Максвелла. Используемый в работе алгоритм состоит в применении комбинации конечно-разностного метода во временной области (FDTD-метод) для решения уравнений Максвелла, метода идеально согласованного слоя — PML-метода для реализации граничных условий и метода общего и рассеянного полей (TF/SF, total field / scattered field) для возбуждения источника поля.
Рассмотрим численный алгоритм решения задачи на примере двумерной прямоугольной области 5 = {(х,г): 0 х ^ а, (рис. 1).
х а *2
TF/SF-трюшци Д_
PML-слои __
Н
□ □ □ □ □ □
□ □ JH □ □ □
£1
□ □ □ □ □ □
□ □ □ □ □ □
□ □ □ □ □ □
□ □ □ □ □ □
-------JVpc-
0 Ъ X
Рис. 1. Схема расчетной области
В основе алгоритма лежит РОТБ-метод [3, 4]. Введем в заданной области разностную сетку
ши = {*/ = ^х ■ г! = А? ■ ¿ = 0,1.....Ых, / = 0,1,...
. .. Мг} [5]. По временной переменной также вводится сетка Шх = {4- = ят, я = 0,1,...}. Значение сеточной функции Т7 в некотором узле сетки (л';,2,) в момент времени Ь обозначим через = Р(1Ьх,]Ьг, ят). В соответствии с РОТО-методом компонента поля Еу рассматривается в узловых точках области (/,/) и в моменты времени а компоненты Нх, Нг — в промежуточных точках области (/,/ + 1/2), (/ + 1/2,/) и в промежуточные моменты времени 5+1/2, т.е. в дробных узлах и на дробных временных слоях. Это приближение позволяет записать пространственные и временные производные в уравнениях Максвелла через центральные разности и получить разностную схему для решения.
Например, для компоненты Еу
(т1 = шь+т-Ьг i WJ' 1/2 - !
S"f" 1 /2 ',/-1/2
c^f'i1/!,/1/2
' —1/2,/ ■
ч)
Для ограничения области в построенном алгоритме используются поглощающие граничные условия, для постановки которых используется метод идеально согласованного слоя — PML-метод (perfectly matched layer) [6, 7]. В этом методе каждая компонента электромагнитного поля разделяется на две части, которые удовлетворяют определенным уравнениям. Например, в двумерном случае ГМ-поля для компоненты Еу
дЕ,, _ дНх
е°е~т~+ ^ =
Щ* , р дНг
Значения параметров ах, аг, а*, а* зависят от расположения точки в области: внутри PML-слоя — области вблизи границ — они отличны от нуля и увеличиваются по мере приближения к границе расчетной области, а в остальной части расчетной области равны нулю.
Для возбуждения волноведущей системы в данном алгоритме используется метод TF/SF [8]. Этот метод состоит в разделении расчетной области на две подобласти: область общего поля и область рассеянного поля. Так как уравнения Максвелла линейны, электромагнитное поле можно представить в виде суммы падающего и рассеянного полей, которые также будут удовлетворять уравнениям Максвелла.
В области рассеянного поля рассматривается только рассеянное поле, а в области общего поля — рассеянное и падающее. На границе между этими областями ставятся специальные граничные условия. Например, для ТЕ -случая эти условия имеют следующий вид для левой границы г^п = Лг/ieft:
(£</,tot)?j,', = {Еу, tot);,/„„ +
Г/гг \s+l/2 / тт \s+l/2 / тт \s+l/2 1
КПХу m)i jnl , i/2~\nx, scat,»,- /ii[1_1/2-*./7A-, Шс^,- х/2 ■
eaeh?
Задавая падающее поле специальным образом, можно рассматривать рассеяние электромагнитных волн на объекте, расположенном в области общего поля. Падающее поле можно получить, например, решая вспомогательную одномерную задачу, в которой возбуждение поля осуществляется током:
Дл\ t) = ехр(0.2тг(/тах - fmln)t - 2) sin(wf),
где ш = 7г(/тах + /т1п), /max, /тт — изменяемые параметры.
Алгоритм также построен для волноведущих систем, основанных на нелинейных фотонных кристаллах с материальными уравнениями:
D = е0 еЕ + col'l-E]2-^ В = ßoßH.
Применение FDTD-метода приводит к разностным уравнениям, которые не позволяют явно выразить зна-
чения компонент поля на следующем временном слое через значения на предыдущем слое. Поэтому сначала вычисляются компоненты векторов О, В и Н, а затем с помощью метода Ньютона строится решение уравнения для компонент вектора Е [9]. Итерационные схемы для определения этих компонент следующие:
(С 1 1 I* 11 _ / п 111* Кпух)и | — КС-ух)-,^
еье;ЛЕухуЛ 1+ 4^® [(£,.s-)fj 11*]' - (A;,)?j 1
\S111*
у
e0eLi + 12г0\'(3) [{ЕуХГ,
\s I 11* _
Чг>1,! I
е0еи(Еуг)^ 1 + 4e0x(3) [ (^)fj 1 -
IP \S I 1 I* 1 1 _ ( p \s I 1
si 1 /
еаеи + 12еох{3> [(АДИТ
где к — порядок итерации. В качестве нулевого приближения берется значение компонент вектора Е в предшествующий момент времени.
Результаты расчетов двумерных волноведущих систем
В первую очередь построенный алгоритм был применен для расчета прохождения электромагнитного импульса через фотонный кристалл с двумерной структурой (см. рис. 1). Значения параметров расчетной области: а = 1.9-10-5 м, Ь = 2.5-10~5 м, Их = /ь = 5- 1СГ8 м, т = 8 • 10"17 с, Я = 3 • 10^7 м. Значения параметров фотонного кристалла: И = 2 ■ 10~7, а = 5 • 10~7, £1 = 1, е-2 = 11.56, Ырс = 30. Толщины РМЬ-слоев составляют 1.9 • 10"6, 2.5 • 10~б по осям х и 2 соответственно. Падающий электромагнитный импульс характеризуется центральной частотой 0.5с/а и полушириной 0.4с/а. Спектр этого импульса имеет форму, изображенную на рис. 2 пунктирной линией.
—А
—В 1 \
/ щ
1 J к.
1.0
в_ 0.8 0.6 0.4 0.2
О
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
/(с/а)
Рис. 2. Модель двумерного фотонного кристалла и его спектральная характеристика. А — падающий импульс, В — спектр пропускания ФК
Как видно из рис. 2, после прохождения импульса через фотонный кристалл, в спектре появляются запрещенные зоны — области, в которых коэффициент пропускания близок к нулю. Изменение диэлектрической проницаемости кристалла приводит к смещению запрещенных зон. Центру первой запрещенной зоны
соответствует частота Д = -1, где г = г1У+г2(1 — у)< и — относительная доля объема, занимаемого средой с диэлектрической проницаемостью . Для используемых параметров эта частота имеет значение /1 =0.305с/а, что подтверждает рис. 2.
Уменьшим величину диэлектрической проницаемости кристалла — пусть £2 = 8. Тогда частота, соответствующая центру первой запрещенной зоны, увеличится и будет равна Д = 0.343с/а. Спектральные характеристики такого кристалла изображены на рис. 3, из которого видно, что запрещенная зона сместилась вправо — в сторону увеличения частоты. Для кристалла с е-2 = 20 частота, соответствующая центру первой запрещенной зоны, равна Д =0.249 с/а — запрещенная зона находится левее, чем у исходного кристалла.
Запрещенные зоны существуют для различных поляризаций. Описанные выше вычисления сделаны для случая ТЕ-поляризации (вектор электрического поля направлен параллельно осям стержней). На рис. 4 изображены также спектральные характеристики фотонного кристалла в ТМ -случае (вектор магнитного поля направлен параллельно осям стержней). Из рисунка видно, что электромагнитные импульсы, имеющие различные поляризации, после прохождения через фотонный кристалл имеют спектры, в которых есть общие
1.5
—А
—В (е2 = 20) —В (е2 = 11.56) ......В(е 2 = 8)
0.5 /(с/а)
Рис. 3. Спектральные характеристики фотонных кристаллов с различными значениями диэлектрической проницаемости
1.0 0.8 0.6 0.4 0.2 0
........А
-В (ТЕ)
----В (ТМ) • \
5\А и
V
I1'!/1 'к1
II1'/ 1 ,1; ! 1 1 1 1
!' I // ' 1 [ V '41 1 [и
/У , ) 1 1 ■ 1 Г\\
1 1 и о ' 1 (-'V А--.......
0.2
0.4
0.6
0.8 Ас/а)
запрещенные зоны. Частоты, для которых запрещенные зоны совпадают, образуют полную запрещенную зону — в этом частотном диапазоне свет любой поляризации не может распространяться в ФК ни в каком направлении.
Удаляя ряд стержней из фотонного кристалла, можно создавать фотонно-кристаллические волноводы. Благодаря наличию в фотонном кристалле запрещенных зон такие волноводы будут работать в диапазоне, совпадающем с диапазоном запрещенных зон. На рис. 5 показаны спектральные характеристики простых вол-новедущих систем, основанных на фотонных кристаллах — волновода, изгиба, разветвителя.
А\
\В
0.35 /(с/а)
Рис. 4. Спектральные характеристики фотонного кристалла для ТЕ - и ТМ-волн
Рис. 5. Спектральные характеристики фотонно-крис-таллического волновода, изгиба и разветвителя
С помощью построенного алгоритма в рамках настоящей работы было проведено моделирование ФК Я-фракталов. Спектральные характеристики фотонных кристаллов Я-фрактала продемонстрированы на рис. 6.
Алгоритм расчетов для трехмерных систем
По аналогии с описанным выше алгоритмом для двухмерной области строится алгоритм для расчета волноведущих систем в трехмерной области. Задача для трехмерной области является достаточно ресурсоемкой, поэтому для ее решения используется параллельное программирование.
В
0.8/(с/а)
б
В
0.8/(с/а)
В
0.8/(с/а)
Рис. 6. Спектральные характеристики фотонного кристалла с обычной структурой (а)
и в виде Я-фракталов (б, в)
Вычисления в работе проводились с помощью вычислительного кластера НИВЦ МГУ на суперкомпьютере СКИФ МГУ «Чебышев». Использовалась наиболее распространенная технология программирования для параллельных компьютеров с распределенной памятью — MPI (message passing interface).
Схема распараллеливания расчета следующая. Трехмерная расчетная область разделяется на несколько одинаковых частей, имеющих общие граничные слои, в соответствии с количеством задействованных в расчете процессоров так, как показано на рис. 7. На каждом временном шаге сначала автономно решается задача в каждой расчетной области на отдельном процессоре, затем полученное решение для точек в граничных слоях передается процессорам, в которых рассчитываются соседние области. Таким образом, расчет электромаг-
нитного поля на каждом шаге по времени состоит из двух этапов: расчета, который происходит параллельно во всех областях, и передачи данных в граничных слоях, которая происходит последовательно.
Такое распараллеливание позволяет уменьшить время расчетов в количество раз, практически равное числу задействованных процессоров (если не учитывать время, затрачиваемое на передачу данных, которое, как показывают расчеты, мало для небольшого количества процессоров). На рис. 8 показана зависимость времени расчетов от количества задействованных процессоров. Эта зависимость соответствует закону Амдала для ускорения:
1 р
%У
Ж
ж
п.
п,
п
N
Рис. 7. Схема разделения расчетной области на несколько областей в зависимости от количества задействованных процессоров N. Серым цветом отмечены общие граничные слои для соседних областей
Рис. 9. Сравнение расчетного и теоретического коэффициента отражения от границы двух различных диэлектрических сред
Ускорение расчета
Время расчета
1 2 3 4 5 Количество процессов
Рис. 8. Зависимость времени расчета и коэффициента ускорения от количества процессоров
где / — доля операций, которые выполняются последовательно, р — число процессоров [10].
В работе расчеты проводились с использованием 20-32 процессоров в зависимости от задачи, расчетная область при этом составляла до 600 х 800 х 600 точек.
Результаты расчетов для трехмерных волноведущих систем
В качестве задачи для тестирования программы рассмотрена задача отражения волны от границы двух различных диэлектрических сред в трехмерной области. На рис. 9 приведена схема расчетной области и график зависимости рассчитанного коэффициента отражения от вставки в сравнении с аналитическим решением. В качестве аналитического решения используется следующая известная зависимость:
к,=
•х/гГ
или, учитывая, что в расчетах = 1,
Сравнение показывает хорошее соответствие полученного результата и теоретического значения.
В трехмерной области рассмотрена задача прохождения электромагнитного импульса через фотонный
Лс/а)
Рис. 10. Сравнение спектров падающего поля (А), рассеянного поля в случае расчетов в двумерной области (В'2) и рассеянного поля в случае трехмерной области (Ш)
1.0
0.8
0.6
0.4
0.2
Рис. И. Спектр падающего поля (А) и рассеянного поля трехмерным фотонным кристаллом (В). На вставке справа изображен пример структуры такого кристалла
0.2 0.4 0.6 0.8 /(с/а)
.......А
-В
кристалл с двумерной структурой (см. рис. 2). Сравнение полученного результата с аналогичным расчетом в двухмерной области показано на рис. 10. Из рисунка видно, что в трехмерной области при прохождении импульса через такую систему рассеяние значительно больше по сравнению с аналогичным расчетом для двумерной области.
На рис. 11 показаны результаты расчетов для трехмерного кристалла. Структура такого кристалла изображена на этом рисунке справа. В расчете использовались следующие параметры: размер стержней, образующих структуру кристалла, составляет 3 • 10^6 х 3 • 10^6 х 90 • 10^6 м, количество слоев — 30, их диэлектрическая проницаемость равна 12.25. Падающее поле направляется на эту структуру сверху.
Заключение
Проведенные в настоящей работе исследования показали, что построенный алгоритм может быть эффективно использован при расчете разнообразных волно-ведущих систем в двумерных и трехмерных областях, в том числе основанных на фотонных кристаллах.
Построенная математическая модель достаточно полно и адекватно описывает физические процессы, происходящие в подобных волноведущих системах. Разработанная на основе данной методики программа
обладает достаточной гибкостью и позволяет вычислять спектральные характеристики устройств, основанных на фотонных кристаллах, с различной структурой. Применение комбинации PML и TF/ST методов дает возможность существенно сократить вычислительные затраты за счет уменьшения расчетной области. Использование распараллеливания расчетов также существенно сокращает время вычислений.
Список литературы
1. Yablonovitch Е. // Phys. Rev. Lett. 1987. 58. P. 2059.
2. Наний O.E., Павлова Е.Г. // Lightwave Russian Edition. 2004. № 3. C. 47.
3. Yee K.S. // IEEE Trans. Antennas Propagat. 1966. 14. P. 302.
4. Taflove A., Hagness S.C. Computational Electrodynamics. Norwood, MA, 2000.
5. Самарский A.A. Теория разностных схем. М., 1983.
6. Ильгамов М.А., Гильманов А.Н. Неотражающие условия на границах расчетной области. М., 2003.
7. Berenger J.-P. // J. Comput. Phys. 1996. 127. P. 363.
8. Anantha V., Taflove A. // IEEE Trans. Antennas Propagat. 2002. 50, N 10. P. 1337.
9. Калиткин H.H. Численные методы. M., 1978.
10. Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. СПб., 2002.
Research of propagation of electromagnetic pulses through photonic crystal structures A.N. Bogolyubov", I.A. Butkarev, Yu.S. Dementieva
Department of Mathematics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia. E-mail: "bogan@yandex.ru.
In this paper we present mathematical model and numerical algorithm allowing to investigate two-dimensional and three-dimensional photonic crystals and waveguiding systems based on photonic crystals. We demonstrate spectral characteristic of photonic crystals and systems based on them.
Keywords: photonic crystal, numerical modeling, parallel computing. PACS: 02.60.Cb. Received 8 July 2010.
English version: Moscow University Physics Bulletin 6(2010).
Сведения об авторах
1. Боголюбов Александр Николаевич — докт. физ.-мат. наук, профессор, профессор; e-mail: bogan@yandex.ru.
2. Буткарев Иван Андреевич — канд. физ.-мат. наук, ст. науч. сотр.; e-mail: butkarev@physics.rnsu.ru.
3. Дементьева Юлия Сергеевна — аспирантка; e-mail: jdementieva@gmail.com.
