Численное моделирование двумерных фотонных кристаллов

В конце 80-х годов прошлого века была впервые предложена концепция фотонных кристаллов (ФК) [1].

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

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

Второе свойство отличает фотонный кристалл от обычной дифракционной решетки. Оно означает, что в определенном спектральном диапазоне свет любой поляризации не может распространяться в ФК ни в каком направлении в одном, двух или трех измерениях. Это и есть уникальное свойство фотонного кристалла, с которым принято связывать возможные революционные события в технике оптической связи, физике лазеров и оптической компьютерной технологии [2].

Фотонные кристаллы предоставляют новые возможности управления световыми потоками благодаря наличию полной запрещенной зоны в плотности электромагнитных состояний в заданной области частот. Внесение искажений в структуру кристалла позволяет создавать различные устройства, например волноводы, резонаторы, фильтры, разветвители, и др. [3,4,5,6,7].

Из-за отсутствия природных фотонных кристаллов первоочередной задачей является создание искусственных фотонно-кристаллических структур [8,9].

В данной работе рассматривается двумерный фотонный кристалл, диэлектрическая проницаемость в нём периодически меняется в двух направлениях. Он представляет собой набор диэлектрических стержней выставленных в узлах прямоугольной решетки с шагом H (рис. 1). Толщина стержней составляет h, а диэлектрическая проницаемость — e₂. Удалив ряд стержней в таком кристалле можно создать фотонно-кристаллический волновод.

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

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

Рассмотрим уравнения Максвелла в ТМ случае (E_x=E_z=H_y=0) в отсутствие токов и свободных зарядов:

Задача решается в прямоугольной двумерной области S={(x, z): 0£x£a, 0£z£b}. Введём в этой области разностную сетку [12] w_h={x_i=ih_x, z_j=jh_z, i=0,1,…,N_x, j=0,1,…,N_z}. Отрезок времени, на котором рассматривается решение задачи, также разбивается на конечное число интервалов — вводится сетка по временной переменной w_s={t_s=st, s=0,1,…}. Значение сеточной функции F в некотором узле сетки (x_i, z_j) в момент времени t_s обозначим через . Компонента поля E_y рассматривается в узловых точках области (i, j)и в моменты времени s, а компоненты H_x, H_z — в промежуточных точках области (i, j+0.5), (i+0.5, j) и в промежуточные моменты времени s+0.5.

Такое рассмотрение позволяет записать пространственные и временные производные в уравнениях через центральные разности. Например, для функции F в точке x_i=ih_x, z_j=jh_z в момент времени t=st производные по x и по t будут иметь вид:

Представляя в таком виде все производные в уравнениях (1)-(3), их можно переписать в следующем виде:

Уравнения (6)-(8) позволяют явно выразить значения компонент электромагнитного поля на следующем временном слое через значения на текущем слое:

Аналогично рассматривается ТЕ случай (H_x=H_z=E_y=0). Уравнения Максвелла в отсутствие токов и свободных зарядов имеют вид

Разностная схема в данном случае принимает следующий вид:

Данные схемы имеют второй порядок аппроксимации по времени и пространственным координатам. Условие сходимости этих схем:

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

Для постановки поглощающих граничных условий используется PML-метод (Perfectly Matched Layer) [14].

В соответствии с данным методом каждая компонента электромагнитного поля разделяется на две части. В декартовых координатах 6 компонент поля E_x, E_y, E_z, H_x, H_y, H_z заменяются на 12 подкомпонент: E_xy, E_xz, E_yx, E_yz, E_zx, E_zy, H_xy, H_xz, H_yx, H_yz, H_zx, H_zy, удовлетворяющих уравнениям:

В двумерном случае, из-за симметрии задачи, получается четыре уравнения и только одна компонента разделяется на две части. В случае TM поля разделяется компонента E_y на подкомпоненты E_yx, E_yz и рассматриваются уравнения

В случае TE поля компонента H_y разделяется на H_yx, H_yz и рассматриваются уравнения

Применяя к уравнениям (12)-(15) и (16)-(19) FDTD-метод, получим соответствующие разностные схемы (20)-(23) и (24)-(27):

Чтобы получить значения исходных компонент поля, нужно сложить значения соответствующих подкомпонент:

Значения параметров зависят от расположения точки в области: внутри PML-слоя они отличны от нуля, а вне — равны нулю (рис. 2).

Рассмотрим двумерную прямоугольную область (рис. 3).

Рис. 3. Схема области, в которой исследуются граничные условия PML

На границах z=0 и z=b введены PML-слои. Значения параметра s_z заданы следующим образом:

Численно моделирование показало, что при увеличении ширины N₀ PML-слоя коэффициент отражения от границы уменьшается. При увеличении параметра s_max до определённого значения коэффициент отражения уменьшается, а при дальнейшем увеличении — растёт. В таблицах 1 и 2 приведены результаты измерений коэффициента отражения от границы двумерной прямоугольной области с размерами a=b=10, N_x=N_z=500, h_x=h_z=0.02 для различных значений параметров PML-слоя s_max и ширины N₀:

s_max	N₀	R	s_max	N₀	R
10	10	0,482	10	50	0,024
10	20	0,205	10²	50	0,007
10	50	0,029	10³	50	0,002
10	100	0,015	10⁴	50	0,011
10	200	0,011	10⁵	50	0,037
Таблица 1. Зависимость коэффициента отражения R от ширины PML-слоя			Таблица 2. Зависимость коэффициента отражения R от параметра s_max

Для возбуждения волноведущей системы используется метод total-field/scattered-field (TF/SF) [15]. Данный метод состоит в разделении расчетной области на две подобласти: область общего поля и область рассеянного поля (рис. 4).

Рис. 4. Схема разделения расчетной области при применении метода TF/SF

Так как уравнения Максвелла линейные, то электромагнитное поле можно представить в виде суммы падающего и рассеянного полей:

где индекс пад означает падающее поле, рас — рассеянное поле, общ — общее поле.

Область общего поля выбирается так, чтобы в ней полностью находился рассеивающий объект (существуют варианты TF/SF метода, в которых преодолено это ограничение), вся оставшаяся часть — это область рассеянного поля. В области рассеянного поля рассматривается только рассеянное поле, а в области общего — рассеянное и падающее. На границе между этими областями ставятся специальные граничные условия. Например, для одномерного TМ-случая эти условия на левой границе между областями общего и рассеянного полей имеют вид:

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

Индексы j_лев, j_прав, i_нижн, i_верхн соответствуют левой, правой, нижней и верхней координатам границы между областями общего и рассеянного полей.

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

В данном случае возбуждение поля осуществляется при помощи заданного тока:

где , , t₀=2pw, w=p(f_max+f_min), f_max, f_min — изменяемые параметры. Спектр импульса, образованного таким током меняется в частотной области от f_min до f_max.

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

Двумерный фотонный кристалл рассматривается в области S={(x, z): 0£x£a, 0£z£b} (рис. 5). В этой области моделируется электромагнитный импульс, распространяющийся в положительном направлении оси z и падающий на кристалл.

Рис. 5. Область, в которой проводились численные расчеты

Кристалл характеризуется параметрами H — постоянная решётки кристалла, ε₁, ε₂ — диэлектрические проницаемости, N_pc характеризует размер кристалла. Во всех измерениях диэлектрические проницаемости ε₁ и ε₂ равны соответственно 1 и 8, h=0.05, H=0.1, N_pc =40. Размеры расчетной области составляют: a=8, b=10. Расчеты велись на разностной сетке c h_x=h_z=0.025, N_x=320, N_z=400 и шаге по времени t=0.012. На границе расчетной области поставлены поглощающие граничные условия типа PML с параметрами: s_max=1000, толщина PML-слоя составляет 10% от размеров области: рядом с границами z=0 и z=b — 1, x=0 и x=a — 0.8. Область общего поля имеет размеры 6.3×6.3 и расположена симметрично в центре расчетной области. Параметры, определяющие спектр импульса, падающего на кристалл, выбраны следующим образом f_min=0, f_max=10.

После прохождения импульса через фотонный кристалл, в спектре появляются запрещённые зоны — области, в котороых коэффициент пропускания близок к нулю. При увеличении размера кристалла N_pc эти зоны имеют более чёткий вид (рис. 6).

Рис. 6. Спектр поля, прошедшего через фотонный кристалл: а — N_pc=10, б — N_pc=50

Запрещённые зоны существуют для различных поляризаций. Электромагнитные импульсы, имеющие различные поляризации, после прохождения через фотонный кристалл имеют спектры, в которых есть общие запрещённые зоны. На рис. 7 показаны соответствующие спектры для TМ-поля (E_x=E_z=H_y=0) — рис. 7а и TE-поля (H_x=H_z=E_y=0) — рис. 7б.

Рис. 7. Спектры поля разной поляризации, прошедшего через фотонный кристалл: а — ТМ, б — ТЕ

Удаляя ряд стержней из фотонного кристалла можно создать фотонно-кристаллический волновод, работающий в диапазоне запрещенных зон (рис. 8).

Рис. 8. Модель фотонно-кристаллического волновода

Из рис. 9 видно, что падающая на границу фотонного кристалла плоская волна частично втягивается в волновод с области большей, чем сечение волновода. Следовательно, в данном случае, управляя характеристиками фотонного кристалла на границе (положением стержней), можно изменять степень проникновения электромагнитного поля в волновод.

Рис. 9. Спектр поля внутри фотонно-кристаллического волновода: а — в точке А, б — в точке В

Управлять направлением распространения света можно при помощи волноводных изгибов.

На рис. 10 изображена модель фотонно-кристаллического волновода с изгибом.

Рис. 10. Модель фотонно-кристаллического волновода с изгибом

Спектр рассеянного поля на выходе из этого изгиба изображён на рис. 11а. Если сравнить этот спектр со спектром поля, рассеянного аналогичным кристаллом без изгиба (рис. 7а), то видно, что в запрещённой зоне возникла полоса пропускания.

Из рис. 7б и рис. 11б видно, что аналогичные свойства наблюдаются и для другой поляризации.

Рис. 11. Спектр поля в изгибе фотонно-кристаллического волновода (точка В): а — ТМ-поле, б — ТЕ-поле

Еще одним вариантом волновода, основанного на фотонных кристаллах, является волновод связанных резонаторов (CCW) рис. 12а. На рис. 12б и 13б приведены вычисленные характеристики этого волновода и изгиба такого волновода, параметры фотонного кристалла и разностной сетки оставлены прежние.

Рис. 12. Результаты расчета CCW: а — схема CCW, б — результаты

Рис. 13. Результаты расчета изгиба CCW: а — схема изгиба, б — результаты

На рис. 14 приведены результаты расчета разветвления волновода связанных резонаторов.

Рис. 14. Результаты расчета разветвления CCW: а — схема разветвления, б — результаты

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

Полученные численные результаты для волноводного изгиба показывают, что фотонные кристаллы являются очень удобными объектами для управления распространением света. Фотонные кристаллы позволяют существенно уменьшить размеры устройств, управляющих распространением электромагнитных волн, и являются перспективными объектами для применения в интегрально оптических устройствах. Рассмотрение задачи во временной области позволяет сразу получить характеристики волноведущей системы в определённом диапазоне частот, определяемом спектром падающих волн.