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

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

А. Н. Боголюбов1, Г. В. Белокопытов2, Ж. О. Домбровская1^

Московский государственный университет имени М. В. Ломоносова, физический факультет, 1 кафедра математики; 2 кафедра физики колебаний. Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2.

E-mail: a dombrovskaya@physics.msu.ru

Статья поступила 23.04.2013, подписана в печать 29.05.2013.

Разработан и реализован численный алгоритм, основанный на совместном применении метода конечных разностей во временной области (FDTD) с поглощающими граничными условиями (комбинация идеально согласованного слоя (PML) и условия Мура первого порядка) и метода полного и рассеянного полей (TF/SF). Приводится полная формулировка алгоритма. Представлены результаты вычислений спектральных характеристик для двумерного фотонного кристалла с дефектами в виде волноводов, которые демонстрируют возможность производить вычисления в широком диапазоне частот.

Ключевые слова: фотонный кристалл, фотонно-кристаллические волноводные системы, запрещенные зоны, метод конечных разностей во временной области (FDTD), идеально согласованный слой (PML), метод разделения полей (TF/SF).

УДК: 519.63. PACS: 02.70.Bf.

Введение

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

Теоретическое исследование спектральных характеристик прошедшего излучения может проводиться аналитически или численно. Среди методов численного моделирования выделяется метод конечных разностей во временной области (FDTD — Finite Difference Time Domain), учитывающий сложную геометрию моделируемого объекта и позволяющий решать нестационарные задачи [2, 3]. Трудности с реализацией данного метода возникают из-за нефизичного переотражения волн от границы расчетной области. В настоящей работе подробно изложен численный алгоритм, основанный на комбинации метода FDTD при наличии идеально согласованного слоя (PML — Perfectly Matched Layer) для учета эффектов у границы расчетной области с методом полного и рассеянного полей (TF/SF —

Total-field/Scattered-field) для моделирования источника в виде плоской волны. Совместное применение PML и TF/SF методов дает возможность сократить вычислительные затраты за счет уменьшения требуемой расчетной области. В настоящей работе приводятся и обсуждаются примеры численных расчетов спектральных зависимостей, получаемых в результате прохождения электромагнитного импульса через двумерный конечный ФК и волноводные системы (прямой волновод и волновод с изгибом) на его основе.

1. Физическая постановка задачи

В общем случае структура ФК (рис. 1) представляет собой пластинку с диэлектрической проницаемостью £pi, в которой имеется N периодов, состоящих из N цилиндрических отверстий радиусом r, заполненных диэлектриком ег. Центры отверстий находятся в узлах квадратной решетки с периодом a в плоскости xy . Слева из внешней неограниченной среды на структуру

N

О

Кр,

Ш

Ф

Рис. 1. Структура двумерного фотонного кристалла

нормально падает квазимонохроматическии импульс с максимальной амплитудой, равноИ единице.

Электромагнитное поле как внутри ФК, так и вне его описывается системой уравнений Максвелла. В работе исследуется прохождение электромагнитной волны с Т£-поляризацией (компонента вектора напряженности электрического поля Е перпендикулярна плоскости распространения волны) через ФК. В декартовой системе координат для однородной немагнитной среды при отсутствии свободных зарядов и токов система уравнений имеет следующий вид:

= c

dt

dHx dt

dHy = c dEz

=c

dHy,

dx dEz

dHx дУ

дУ '

dt dx

где c — скорость света в вакууме. В случае линейной, изотропной и бездисперсионной среды вектор электрического смещения D связан с вектором напряженности электрического поля E материальным уравнением

Dz = e(x, y)Ez, (2)

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

2. Алгоритм численного моделирования 2.1. Метод конечных разностей во временной области

Для численного решения системы уравнений (1) используется метод FDTD [4]. Он является одним из наиболее популярных методов численного решения нестационарных задач для системы вихревых дифференциальных уравнений Максвелла. Идея этого метода основывается на алгоритме Йе [5], согласно которому сетки для электрического и магнитного полей (рис. 2) выбираются так, чтобы они были равномерны и смещены по отношению друг к другу на половину шага дискретизации по каждой из пространственных переменных и по времени. Введенная таким образом сетка наиболее адекватно соответствует физике процесса: изменение электрического поля E во времени, как и в уравнениях Максвелла, зависит от его значения на предыдущем шаге и от пространственного изменения

J+ 1/2 -----Ж"-----

Я„

Яг

—-ж-—

У

j-1/2 -

-ж---

! * i - 1/2

Яг

нг

---ж-

i + 1/2

2 i

Рис. 2. Сетки для электрического и магнитного полей

магнитного поля Н, аналогично изменение поля Н зависит от распределения в пространстве поля Е.

Для построения консервативной разностной схемы используется интегро-интерполяционный метод [6], вводятся дробные узлы (на рис. 2 отмечены крестиками). При таком рассмотрении приближенные пространственные и временные производные в системе (1) можно записать через центральные разности. Для производных по координате (например, по х) и времени в точке (г,]) в п-м шаге по времени имеем

du n dx

uli+1/2j u U—1/2,/ + -" + 1/2,/ " H —1/2,j

, ч du

(1) ~m

i,/ hx

n u|n+ 1/2 - u|n— 1/2

i,/

T

+ O(T2)

u|n+1/2,/ u|i—1/2,/ ~ hx

u|n+1/2 - u|n—1/2

(3)

где Нх и т — шаги сетки. Используя формулы типа (3), перепишем систему уравнений (1) в конечно-разностном виде

ог1п+1 = +

+ CT

H in+1/2 _ H in+1/2 H in+1/2 _ H in+1/2

Hy |i+1/2,/ Hy |i—1/2,/ Hx|i,/+1/2 Hx|i,/—1/2

hx

Hx|n:/++1/22=Hx - C

E n E n Ez|i,/+1 EzU,i

(4)

H in+1/2 = H in—1/2 ПУ | i+1/2,/ = ПУ U+1/2/

+ CT

E n E n Ez|i+l,/ EzU

i,j

hx

Построенная схема является явной и, как хорошо известно, обладает условной устойчивостью. Ограничение на выбор значений шагов по пространству и времени накладывается условием Куранта [7]. В двумерном случае оно имеет следующий вид:

< 1.

h

(5)

Здесь и в дальнейшем рассматривается квадратная сетка по координатам (Нх = Ну = Н). Отметим, что в случае равенства в (5) схема становится неустойчивой из-за анизотропии численной дисперсии. Анизотропия проявляется как результат того, что в двумерной задаче излучение рассеивается под различными углами, и рассеянным волнам приходится преодолевать разные расстояния для достижения границы расчетной области. Исследование [8] показало, что максимальное значение временного шага, при котором система (4) устойчива, составляет

т- = 6 СН/2 " 06 Н (6)

Для простоты вычислений и наглядности получаемых результатов коэффициент в (6) удобно выбрать равным 0.5.

2.2. Метод разделения полей

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

y

h

y

функцию с несущей частотой /0 = /тах+ т оптического на верхней границе (/ = /ь) сигнала, модулированную гауссовым профилем в направлении распространения. Такой импульс задается

Г) |"+1 _ Г) п ст тг¡пе|"+1/2

Б:'г,/ь = Б:'г,/ь - |у£) +1/2,

упрощением системы (1) до одномерного случая:

дБ:

дНХ

дt

дНХпс

дt

= с-

= — с-

дх дЕ1пе дх

+ /(0,

Н |"+1/2 = н I"—1/2

Нх |г,/ь + 1/2 = Нх1;,/6 + 1/2

(7) на левой границе (г = га)

причем возбуждение поля осуществляется током в виде функции Габора

¡(t) = si n(2п/ot) ехр< -

И ■

(8)

где

^ = 2а, а = ^Лт (/тах - /тт) , /тах и /тш - изменяемые параметры.

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

Одним из широко распространенных вариантов моделирования плоской электромагнитной волны является использование метода разделения вычислительной области на две подобласти (рис. 3) - полного и рассеянного полей (ТР/БР) [9]. Полное поле состоит из суммы поля падающей волны и рассеянного излучения = Епе + £^еа', Я404 = НШе + Нsеat), а рассеянное поле — только из рассеянного излучения. На границах между ними ставятся специальные условия: на нижней границе (/ = /а)

ст,

1пе|"+1/2

ja-1/2,

Б \п+1 = Б I". +— н

Н "+1/2 = н \п~ 1/2 + ст Е'пе I"

Нх 1 ¡,1а - 1/2 = Нх 1 ¡,1а - 1/2 + кЕ: '¡а ,

Миг

ХР

РМЬ

«а Ч

Рис. 3. Схематичное изображение расчетной области

ст

_ Еше I"

ЙЕ: '¡ь,

ст

Н |"+1/2 = н I"—1/2 _ ^Е 1пе"

НУ '¡а 1/2,/' = ^ '¡а —1/2,/- -

на правой границе (¡ = ¡ь)

Н |"+1/2 = н п—1/2

ПУ Чь + 1/2,/ = ПУ Чь + 1/2,/ ' ^

ст

+ - ЕГ|".

Иными словами, за пределами границы из полного поля вычитается поле падающей волны и для каждого вектора учитывается, какое направление выбрано за положительное. Существенным моментом является то, что падающая волна не взаимодействует с внешней границей 7 расчетной области, на которой ставятся дополнительные условия. Таким образом, метод ТР/БР позволяет создать аналогичные эксперименту условия при численном моделировании процесса.

2.3. Поглощающие граничные условия

При использовании метода РЭТЭ возникает проблема ограничения области, в которой строится решение. Одним из способов численного моделирования ухода волны на бесконечность является использование поглощающих граничных условий.

Вводится понятие идеально согласованного слоя (РМЬ), который является слоем из гипотетического материала [10]. Он представляет собой поглощающую приграничную область, состоящую из некоторого числа слоев, и позволяет получить низкие коэффициенты отражения от границы 7. Работа классического РМЬ основывается на разбиении компонент исходных полей Е и Н на две, для каждой из которых решаются свои уравнения, составляющие систему уравнений Беренге-ра [11].

Рассмотрим распространение в среде 1 плоской волны, падающей на некоторую среду 2, граница раздела между средами выбирается плоская. Величины коэффициентов отражения й = и пропускания

Т = выражаются через волновые импедансы пер-

вой п = л/ 1т и второй среды ц2 = ,/ее* соответ-

у & 1 у е2е

ственно, где е* и ц* — фиктивные диэлектрическая и магнитная проницаемости [12], которые являются комплексными величинами. Вернемся к системе уравнений

Максвелла и, переходя в уравнениях (1), (2) в частотна

ную область, добавим электрические е** = ер + -ш и магнитные = 1 + ^н потери поглощающих слоев:

шОге*(х)е*(у) = с

дну днх дх

дУ

- дЕ:

1Шнх ^*х(х)^х(у) = -с -у,

- дЕ: шну^*у(х)^*у (У) = с -¡х,

Бг(ш) = е(ш)Ег (ш),

где ш — круговая частота; положительные величины а^ и аН — удельные проводимости, определяющие скорость затухания электрической и магнитных компонент соответственно, р = х или у; 1 — мнимая единица.

Существует два условия, которые, как показано в работе [13], необходимо наложить, чтобы сформировать РМЦ полностью поглощающий все падающие на него волны независимо от угла падения и длины волны: 1) импеданс, попадающий на поглощающие слои из внешней среды должен быть постоянной величиной (в данном случае нормирован на единицу)

/до* .

П0 =Пр =у ^ =1;

2) фиктивные диэлектрические проницаемости являются взаимно обратными величинами в том смысле, что для них выполняется следующее условие:

е* = —

х е*'

еу

Аналогичное утверждение имеет место и для фиктивных магнитных проницаемостей:

1

А* =

А

У

С учетом перечисленных условий можно положить ер = 1, а°р = аН = а. (10)

После подстановки (10) в (9) получаем систему

iuDz

шИх 1 +

a(x)

i шИУ

Л + Л + «Ш

У 1ш J \ 1ш J )"('

)('+% )'

= c

H

dx

H dy

ш

1 + ^

i ш

dEz

=c

dy '

1 + ^

=c

i ш

dEz dx

г - (' +

dx \ ш J

-1

dx'

Учет РМЬ сопровождается появлением в системе уравнений (4) дополнительных множителей, зависящих от ширины Ь поглощающей области, обратно пропорциональной величине коэффициента отражения от границы 7 [15], и ограниченных некоторыми величинами с целью уменьшения отражения внутри РМЬ. Положим

= 0.33 ■

L - q + 1 L

е [0,0.33],

где ц = г, ] = 1,..., Ь+1, множитель 0.33 получен эмпирическим путем и обеспечивает стабильную работу разностной схемы при выборе значения величины Ь равным 0.1 от размеров сторон расчетной области,

g^2|i = , 1 , , /¿2|г+1/2 = т+п-е [0.75, 1],

1 + /1|ï

1 + /1|*

+ 1/2

&3\, = , /¿3 \г + 1/2 = е [0.5,1].

1 + /1 \г 1 + /1 \г + 1/2

Аналогичные выражения могут быть записаны для

gj2\j, /¡2\/+1/2, gjз\j, Ш/+1/2. Таким образом, с учетом

РМЬ можно вместо (4) записать

Dz|j1 = giéghljDzlb + g*21 ï gj21 j X

CT

h

YИ in+1/2 _ И n+1/2 \ _ (И in+1/^ _ И n+1/2 H lï+1/2,j ИУ 1 ï —1/2,j^ l4Hx|ï,j+1/2 Hx 1 ï,j — 1/:

Hx|i!++/22 = /i3|ï+1/2Иx |П—+11//2 - /ï2 1 ï +1/2 y [Ez|n,j+1 - Ez ^ ]

ИУ |"-+1/2î,j = /j3 |j+1/2Hy |n+11//2,j + /j3 | j+1/2 h [Ez|n+1,j - ^j] .

Таким образом, РМЬ можно математически трактовать как растяжение или сжатие координат с использованием комплекснозначных функций [14], например:

основанное на аналитическом продолжении уравнений Максвелла в комплексную плоскость таким образом, что их решение экспоненциально затухает.

Для структур, сформированных только из диэлектрических рассеивателей, полного поглощения не наблюдается: возбуждаются артефактные волны, амплитуды которых на несколько порядков меньше чем у падающей волны, но существенны и при отсутствие поглощения могут привести к расходимости РОТО. «Профиль потерь» складывается из отражений от первого слоя РМЬ и между слоями из-за погрешности конечно-разностной дискретизации, а также из отражения от последнего слоя. На границе расчетной области 7 ставятся однородные граничные условия Дирихле:

Ег\7 = 0.

У U+1/2,j J 'У H—1/2,j I V x|ï,j+1/2 J Jx |ï,j—1/2

cT h

+ cT

'У |ï+ 1/2,j / j3 |j+1/2Hy Vï+Hij + /j3 |j+1/2 h

Отметим, что в случае точечного источника достаточно использовать только PML, но в случае разделения полей для задания плоской волны, необходимо наложить дополнительное условие для поглощения падающего излучения в области полного поля после его прохождения через ФК. Причиной этого является расположение PML в области рассеянного поля, где он не может поглотить падающее излучение. Однако учитывая тот факт, что плоская волна, падающая нормально, является по своей сути одномерной (7), в настоящей работе для ее поглощения на поле падающего излучения накладывается условие Мура (Mur condition) первого порядка [16]

Ez|n+1 = Ez|" + CT+h №Î -Ez|?) .

Оно проще для реализации по сравнению с PML и благодаря меньшей затратности позволяет экономить вычислительные ресурсы.

2.4. Получение спектра

С помощью преобразования Фурье после получения значений напряженности электрического поля в некоторой точке области за время счета T можно сразу определить спектральные свойства системы в частотном интервале [/min, /max]

Z(/) = Y,E(t)exp{-2m/T\ .

t=1 ^ '

3. Результаты моделирования

Описанный выше численный алгоритм реализован в виде программы для моделирования прохождения

3

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

соответствующей зоны Бриллюэна [17]:

&

/& * 9 / П )+ , (11)

2 у ер1 (1 - V) + V

где частоты /& измеряются в относительных единицах (с/а), V — относительная доля площади, занимаемая отверстиями, заполненными средой с диэлектрической проницаемостью ег. Формула (11) справедлива в приближении почти свободных фотонов, когда разность |ер1 — ег| мала, для бесконечного ФК в случае предельно длинных электромагнитных волн. Тем не менее она дает удовлетворительное описание и для конечных структур, что позволяет использовать ее с целью предварительного анализа спектральных зависимостей.

Вычисления проводились для ФК, состоящего из 11 х 11 отверстий радиуса г = 0.2а, заполненных диоксидом кремния БЮ2, у которого диэлектрическая проницаемость ег = (1.53)2, и приближение почти свободных фотонов можно считать выполненным. Интерес представляла область видимого излучения (длина волны изменялась в диапазоне от 0.4 до 0.75 мкм). Для заданных параметров центральная частота первой запрещенной зоны равна /1 * 0.463, что достаточно хорошо согласуется с численным экспериментом (рис. 4), в котором получено, что диапазон частот [0.454,0.459] является стоп-зоной в спектре волны с ТЕ-поляри-зацией.

Модифицируя структуру ФК путем удаления части структурных элементов, можно создавать оптические устройства для управления светом, например прямой волновод и волновод с изгибом под углом 90° (рис. 5). Волноводные системы на основе ФК будут работать в частотном диапазоне, совпадающем с областями запрещенных зон у такого кристалла. В случае длительного

f,c/a

Рис. 4. Частотные спектры нормированной интенсивности в точке (¡ь+а, (¡а +/ь)/2) в трех случаях: 1 — для свободного пространства, 2 — для ФК, 3 — для прямого волновода

падающего импульса с узким спектром, локализованным около частоты, соответствующей центру запрещенной зоны, свет практически без потерь распространяется по волноводному каналу. Интенсивность поля на выходе из прямого волновода составляет 0.84 отн. ед., а на выходе из изогнутого волновода равняется 0.79, т. е. потери в изгибе малы и составляют приблизительно 6%. Таким образом, для волноводных систем в области запрещенных частот возникает полоса пропускания, которая в данном случае покрывает всю запрещенную зону (рис. 4). Пунктирной линией изображена зависимость интенсивности падающего импульса, задаваемого в виде (10), от частоты; полужирной и тонкой линиями показаны спектральные зависимости для ФК и прямого волновода, нормированные на единицу. Спектр излучения, прошедшего через волновод с изгибом, не приводится, так как он практически не отличается от спектра для прямого волновода. Результаты вычислений согласуются с работами других авторов, например [18], что свидетельствует о корректности метода.

Заключение

В настоящей работе построена математическая модель процесса прохождения электромагнитного импуль-

ооооооооооо ооооооооооо ооооооооооо • 00.090 000 00

Щ M

Л Л Л Л^^^ г

о о о о о

0 0 0 0 0

о о о о %

о о о о о

о о о о о

10 о о о о

о о о о о

о о о о о

о о о о о

о о о о о

L.0 ).8 ).6 ).4 ).2 )

-0.2 -0.4 -0.6 -0.8 -1.0

о о о о о о

0 0 0 0 0 0

о о о о о о

о о о о о о

о о о о о о

о о о о о о

о о о о о о

о о о о о о

о о о о о о

о о о о о о

о о о о о

0 0 0 0 0

о о о о о

о о о о о

о о о о о

0 0 0 0 0

о о о о о

о о о о о

о о о о о

0 0 0 0 0

Рис. 5. Прохождение квазимонохроматического электромагнитного импульса от точечного источника через фотонно-кристаллические волноводные системы: прямой волновод (слева) и волновод с изгибом под

углом 90о (справа)

са через двумерные фотонно-кристаллические волно-водные системы в широком диапазоне частот. Дано подробное описание численного алгоритма, основанного на комбинации метода FDTD с использованием PML и метода TF/SF. Данный алгоритм позволяет экономить время счета и вычислительные ресурсы. Его работоспособность продемонстрирована на примере расчета спектральных характеристик излучения для ФК на основе диоксида кремния и фотонно-кристаллических волноводных систем.

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

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

1. Yablonovïtch E. // Phys. Rev. Lett. 1987. 58. P. 2059.

2. QïangZ., Zhou W., Sore/ R.A. // Opt. Express. 2007. 15, N 4. P. 1823.

3. Боголюбов А.Н., Буткарев И.А., Дементьева Ю.С. // Вестн. Моск. ун-та. Физ. Астрон. 2008. № 3. С. 33.

4. Sullïvan D.M. Electromagnetic Simulation Using the FDTD Method. N.Y., 2000.

5. Yee K.S. // IEEE Trans. Antennas Propagat. 1966. 14. P. 302.

6. Калиткин Н.Н. Численные методы. СПб., 2011.

7. Kunz K.S., Luebbers R.J. The Finite Difference Time Domain Method for Electromagnetics. Boca Raton (FL), 1993.

8. Fang J. Time-Domain Finite Difference Computations for Maxwell's Equations: Ph.D. diss. Berkeley (CA), 1989.

9. Taflove A., Hagness S.C. Computational Electrodynamics: the Finite Difference Time-Domain Method. Norwood (MA), 2000.

10. Teixeira F.L. // Radio Science. 2003. 38, N 2. P. 8014.

11. Berenger J.P. // J. Comput. Phys. 1996. 127. P. 363.

12. Sullivan D.M. // IEEE Microwave and Guided Wave Letters. 1996. 6, N 2. P. 97.

13. Sacks Z.S., Kingsland D.M., Lee R., Lee J.F. // IEEE Trans. Antennas and Propagat. 1995. 43, N 12. P. 1460.

14. Loh P., Oskooi A.F., Ibanescu M. et al. // Phys. Rev. E. 2009. 79. 065601.

15. Боголюбов А.Н., Буткарев И.А., Дементьева Ю.С. // Журн. радиоэлектроники. 2006. № 5 (http://jre.cplire.ru/jre/nov06/3/text.html).

16. Mur G. // IEEE Trans. Electromagn. Compat. 1981. 23. P. 377.

17. John S. // Phys. Rev. Lett. 1987. 58. P. 2486.

18. Johnson S.G., Joannopoulos J.D. // Acta Materialia. 2003. N 51. P. 5823.

Modelling of spectral dependencies for 2D photonic crystal waveguide systems

A.N. Bogolubov1, G. V. Belokopytov2, Z.O. Dombrovskaya1,a

1 Department of Mathematics; 2 Department of Oscillation Physics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia. E-mail: a dombrovskaya@physics.msu.ru.

The work is devoted to development and realization of the numerical algorithm, based on joint application of the Finite Difference Time Domain method (FDTD) with absorbing boundary conditions (a combination of the Perfectly Matched Layer (PML) and the first order Mur condition) and the Total-field/Scattered-field method (TF/SF). A complete formulation of the algorithm is given. To demonstrate a possibility to make computations in broad frequency range, the results of spectral characteristics' calculation for 2D photonic crystal with waveguide defects are presented.

Keywords: photonic crystal, photonic crystal waveguide systems, band gaps, Finite Difference Time Domain method (FDTD), Perfectly Matched Layer (PML), Total-field/Scattered-field (TF/SF). PACS: 02.70.Bf. Received 23 April 2013.

English version: Moscow University Physics Bulletin 5(2013).

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

1. Боголюбов Александр Николаевич — доктор физ.-мат. наук, профессор; тел.: (495) 939-10-33; e-mail: bogan7@yandex.ru.

2. Белокопытов Геннадий Васильевич — канд. физ.-мат. наук, доцент; тел.: (495) 939-32-61; e-mail: gvb@phys.msu.ru.

3. Домбровская Жанна Олеговна — студентка; тел.: (495) 939-10-33; e-mail: dombrovskaya@physics.msu.ru.