Лицензия MIT
© Алексей Александрович Щербаков, 2024
Модальные методы можно применять для решения различных задач с различной геометрией. Основная идея заключается в том, чтобы рассчитать полный модальный базис в различных областях пространства, а затем произвести сшивку базисов на границах соседних областей при помощи условий сопряжения. При этом, ввиду взаимной неортогональности базисов в разных подпространствах, возникают интегралы перекрытия различных мод. Здесь эта идея будет продемонстрирована на примере слоя одномерного фотонного кристалла.
Воспользуемся результатами модального разложения для одномерного бесконечного фотонного кристалла. Пусть теперь ось $X$ декартовых координат задает направление периодичности. Рассмотрим плоскую границу $z=0$ такую, что полупространство ниже границы при $z<0$ заполнено фотонным кристаллом, а полупространство $z>0$ заполнено однородной изотропной средой с параметрами $\varepsilon$, $\mu$. Зафиксируем блоховский волновой вектор $k_B \equiv k_{x0}$. Разложение поля в верхнем полупространстве по плоским волнам имеет вид \begin{align} F\left(x,z>0,k_{B}\right) &= \sum_{m}\left(a_{m}^{+}e^{ik_{zm}z}+a_{m}^{-}e^{-ik_{zm}z}\right)e^{ik_{xm}x} \\ G\left(x,z>0,k_{B}\right) = \left(\mp\right)_{s,p}\dfrac{1}{i\omega\eta}\left(\dfrac{\partial F}{\partial z}\right) &= \left(\mp\right)_{s,p} \sum_{m}\dfrac{k_{zm}}{\omega\eta}\left(a_{m}^{+}e^{ik_{zm}z}-a_{m}^{-}e^{-ik_{zm}z}\right)e^{ik_{xm}x} \end{align} где $k_{xm} = k_{x0} + 2\pi m/\Lambda$, $k_{zm} = \sqrt{\omega^2\varepsilon\mu - k_{xm}^2}$. Счетный набор плоских волн в разложении возникает вследствие заданной фотонным кристаллом периодичности всего пространства. В области $z<0$ воспользуемся полученным модальным разложением: \begin{align} F\left(x,z<0,k_{B}\right) &= \sum_{m}\left(c_{m}^{+}e^{i\beta_{m}z}+c_{m}^{-}e^{-i\beta_{m}z}\right)v_{m}\left(x\right) \\ G\left(x,z<0,k_{B}\right) &= \left(\mp\right)_{s,p}\sum_{m}\dfrac{\beta_{m}}{\omega\eta\left(x\right)}\left(c_{m}^{+}e^{i\beta_{m}z}-c_{m}^{-}e^{-i\beta_{m}z}\right)v_{m}\left(x\right) \end{align} Непрерывность тангенциальных компонент полей на границе $z=0$ даёт: \begin{align} \sum_{m}\left(c_{m}^{+}+c_{m}^{-}\right) v_{m}\left(x\right) &= \sum_{n}\left(a_{n}^{+}+a_{n}^{-}\right)e^{ik_{xn}x} \\ \sum_{m}\left(c_{m}^{+}-c_{m}^{-}\right)\dfrac{\beta_{m}}{\eta\left(x\right)} v_{m}\left(x\right) &= \sum_{n}\left(a_{n}^{+}-a_{n}^{-}\right)\dfrac{k_{zn}}{\eta_{b}}e^{ik_{xn}x} \end{align} Эту систему можно разрешить, например, для коэффициентов разложения поля в однородной среде. Воспользуемся ортогональностью плоских волн - умножая обе части уравнений на $\exp(-ik_{xp}x)$ и интегрируя по периоду, получаем \begin{align} a_{n}^{+}+a_{n}^{-} &= \sum_{m}\left(c_{m}^{+}+c_{m}^{-}\right)\left[\dfrac{1}{\Lambda}\intop_{0}^{\Lambda}v_{m}\left(x\right)e^{-ik_{xn}x}dx\right] = \sum_{m}\mathcal{I}_{nm}^{(1)}\left(c_{m}^{+}+c_{m}^{-}\right) \\ \dfrac{1}{\eta_{b}}k_{zn}\left(a_{n}^{+}-a_{n}^{-}\right) &= \sum_{m}\beta_{m}\left(c_{m}^{+}-c_{m}^{-}\right)\left[\dfrac{1}{\Lambda}\intop_{0}^{\Lambda}\dfrac{v_{m}\left(x\right)}{\eta\left(x\right)}e^{-ik_{xn}x}dx\right] = \sum_{m}\mathcal{I}_{nm}^{(2)}\beta_{m}\left(c_{m}^{+}-c_{m}^{-}\right) \end{align} где введены обозначения интегралов перекрытия $\mathcal{I}_{nm}^{(1,2)}$. Запишем связь между наборами коэффициентов разложения полей по обе стороны границы раздела полупространств в форме Т-матрицы: \begin{equation} \left(\begin{array}{c} a_{n}^{+}\\ a_{n}^{-} \end{array}\right)=\sum_{m}\left(\begin{array}{cc} T_{11,nm} & T_{12,nm}\\ T_{21,nm} & T_{22,nm} \end{array}\right)\left(\begin{array}{c} c_{n}^{+}\\ c_{n}^{-} \end{array}\right) = \dfrac{1}{2}\left(\begin{array}{cc} \mathcal{I}_{nm}^{(1)}+\dfrac{\beta_{m}}{k_{zn}}\mathcal{I}_{nm}^{(2)} & \mathcal{I}_{nm}^{(1)}-\dfrac{\beta_{m}}{k_{zn}}\mathcal{I}_{nm}^{(2)}\\ \mathcal{I}_{nm}^{(1)}-\dfrac{\beta_{m}}{k_{zn}}\mathcal{I}_{nm}^{(2)} & \mathcal{I}_{nm}^{(1)}+\dfrac{\beta_{m}}{k_{zn}}\mathcal{I}_{nm}^{(2)} \end{array}\right) \left(\begin{array}{c} c_{n}^{+}\\ c_{n}^{-} \end{array}\right) \end{equation} В случае простой ячейки фотонного кристалла, состоящей из двух однородных слоёв, интегралы перекрытия легко находятся аналитически: \begin{align} \mathcal{I}_{nm}^{(1)} =& \dfrac{d_{1}}{\Lambda}\left(\mathcal{A}_{1m}^{+}e^{-\frac{1}{2}i\left(k_{nx}-\kappa_{1m}\right)d_{1}}\mathrm{sinc}\left[\frac{1}{2}\left(k_{nx}-\kappa_{1m}\right)d_{1}\right]+\mathcal{A}_{1m}^{-}e^{-\frac{1}{2}i\left(k_{nx}+\kappa_{1m}\right)d_{1}}\mathrm{sinc}\left[\frac{1}{2}\left(k_{nx}+\kappa_{1m}\right)d_{1}\right]\right) \\ &+ \dfrac{d_{2}}{\Lambda}e^{-ik_{nx}d_{1}}\left(\mathcal{A}_{2m}^{+}e^{-\frac{1}{2}i\left(k_{nx}-\kappa_{2m}\right)d_{2}}\mathrm{sinc}\left[\frac{1}{2}\left(k_{nx}-\kappa_{2m}\right)d_{2}\right]+\mathcal{A}_{2m}^{-}e^{-\frac{1}{2}i\left(k_{nx}+\kappa_{2m}\right)d_{2}}\mathrm{sinc}\left[\frac{1}{2}\left(k_{nx}+\kappa_{2m}\right)d_{2}\right]\right) \end{align}
\begin{align} \mathcal{I}_{nm}^{(2)} = & \dfrac{\eta}{\eta_{1}}\dfrac{d_{1}}{\Lambda}\left(\mathcal{A}_{1m}^{+}e^{-\frac{1}{2}i\left(k_{nx}-\kappa_{1m}\right)d_{1}}\mathrm{sinc}\left[\frac{1}{2}\left(k_{nx}-\kappa_{1m}\right)d_{1}\right]+\mathcal{A}_{1m}^{-}e^{-\frac{1}{2}i\left(k_{nx}+\kappa_{1m}\right)d_{1}}\mathrm{sinc}\left[\frac{1}{2}\left(k_{nx}+\kappa_{1m}\right)d_{1}\right]\right) \\ &+ \dfrac{\eta}{\eta_{2}}\dfrac{d_{2}}{\Lambda}e^{-ik_{nx}d_{1}}\left(\mathcal{A}_{2m}^{+}e^{-\frac{1}{2}i\left(k_{nx}-\kappa_{2m}\right)d_{2}}\mathrm{sinc}\left[\frac{1}{2}\left(k_{nx}-\kappa_{2m}\right)d_{2}\right]+\mathcal{A}_{2m}^{-}e^{-\frac{1}{2}i\left(k_{nx}+\kappa_{2m}\right)d_{2}}\mathrm{sinc}\left[\frac{1}{2}\left(k_{nx}+\kappa_{2m}\right)d_{2}\right]\right) \end{align}Далее перейдем к случаю, когда фотонный кристалл занимает область пространства в пределах плоского слоя толщиной $h$ - $-h/2\leq z\leq h/2$. Обозначим векторы амплитуд разложения поля по плоским волнам в однородных изотропных полупространствах $z<-h/2$ и $z>h/2$ как $\boldsymbol{b}^{\pm}=\left\{ b_{m}^{\pm}\right\} _{m=-\infty}^{\infty}$ и $\boldsymbol{a}^{\pm}=\left\{ a_{m}^{\pm}\right\} _{m=-\infty}^{\infty}$ соответственно, а материальные параметры соответствующих сред как $\varepsilon_{a,b}$ и $\mu_{a,b}$. Амплитуды мод, по которым ведется разложение поля внутри слоя введем следующим образом. Вектор $\boldsymbol{c}^{+}$ будет обозначать амплитуды у нижней границы $z=-h/2+0$, а вектор $\boldsymbol{c}^{-}$ - амплитуды у верхней границы $z=h/2-0$. Считая эти амплитуды самосогласованными, их перенос в пределах слоя осуществляется умножением на вектор фазовых множителей, из которых можно составить диагональную матрицу $\mathcal{E}(\Delta z)=\mathrm{diag}\left\{ \exp\left(i\beta_{m}|\Delta z|\right)\right\} _{m=-\infty}^{\infty}$, где $\beta_{m}$ - константы распространения мод, являющиеся решениями дсперсионного уравнения в бесконечном фотонном кристалле.
Воспользуемся выведенными Т-матрицами границ радела и запишем их для верхней и нижней границ рассматриваемого слоя: \begin{equation} \left(\begin{array}{c} \boldsymbol{b}^{-}\\ \boldsymbol{b}^{+} \end{array}\right) = \left(\begin{array}{cc} T_{11}^{(1)} & T_{12}^{(1)} \\ T_{21}^{(1)} & T_{22}^{(1)} \end{array}\right) \left(\begin{array}{c} \mathcal{E}\boldsymbol{c}^{-} \\ \boldsymbol{c}^{+} \end{array}\right) \end{equation}
\begin{equation} \left(\begin{array}{c} \boldsymbol{a}^{+} \\ \boldsymbol{a}^{-} \end{array}\right)=\left(\begin{array}{cc} T_{11}^{(2)} & T_{12}^{(2)} \\ T_{21}^{(2)} & T_{22}^{(2)} \end{array}\right)\left(\begin{array}{c} \mathcal{E}\boldsymbol{c}^{+}\\ \boldsymbol{c}^{-} \end{array}\right) \end{equation}Эти соотношения можно переписать в форме S-матрицы: \begin{equation} \left(\begin{array}{c} \boldsymbol{b}^{-}\\ \boldsymbol{a}^{+} \end{array}\right) = S \left(\begin{array}{c} \boldsymbol{b}^{+}\\ \boldsymbol{a}^{-} \end{array}\right) \end{equation}
где в явном виде \begin{equation} S = \left(\begin{array}{cc} T_{11}^{(1)}\mathcal{E} & T_{12}^{(1)}\\ T_{12}^{(2)} & T_{11}^{(2)}\mathcal{E} \end{array}\right)\left(\begin{array}{cc} T_{21}^{(1)}\mathcal{E} & T_{22}^{(1)}\\ T_{22}^{(2)} & T_{21}^{(2)}\mathcal{E} \end{array}\right)^{-1} \end{equation} Использование более простой для аналитических вычислений T-матрицы в данном случае приводит к вычислительной неустойчивости и расходимости метода, посколько в матрице переноса появляются множители $\exp\left(-i\beta_{m}|\Delta z|\right)$, которые могут оказаться очень большими для сильно затухающих волн, что приводит к быстрому накоплению ошибки численных расчетов. Затухающие волны же необходимо учитывать при использовании данного численного метода расчета дифракции.
Основными недостатками модального метода, описанного выше, является, во-первых, необходимость решения трансцендентного дисперсионного уравнения, при этом не пропуская его последовательные корни, а, во-вторых, сложность масштабирования подхода на двумерные фотонные кристаллы. Оба недостатка устраняются в рамках так Фурье-модального метода. Последовательность шагов остаётся той же, что при формулировке модального метода, а именно, решение задачи на собственные значения в бесконечном фотонном кристалле, нахождение Т-матрицы плоской границы, и нахождение матрицы рассеяния плоского слоя кристалла. Третий шаг идентичен уже описанному выше, поэтому здесь приведем детали первых двух шагов.
Для всех компонент полей кроме теоремы Блоха запишем разложение в ряд Фурье периодической части волновой функции: \begin{equation} F\left(x,z\right) = \exp\left(ik_Bx\right) F_{\Lambda}\left(x,z\right) = \sum_{m}F_{m}\left(z\right)\exp\left(ik_{xm}x\right) \end{equation} Подставляя подобные разложения в уравнения Максвелла для двумерной задачи и ТЕ-поляризованного поля, получим \begin{align} -\frac{d\boldsymbol{E}_{y}}{dz} &= i\omega\mu_{0}\boldsymbol{H}_{x} \\ K\boldsymbol{E}_{y} &= \omega\mu_{0}\boldsymbol{H}_{z} \\ \dfrac{d\boldsymbol{H}_{x}}{dz} - iK\boldsymbol{H}_{z} &= -i\omega [[\varepsilon]] \boldsymbol{E}_{y} \end{align} где векторы $\boldsymbol{E}_{y}=\left\{ E_{ym}\right\}$, $\boldsymbol{H}_{x}=\left\{ H_{xm}\right\}$, $\boldsymbol{H}_{z}=\left\{ H_{zm}\right\}$, и матрица $K=\mathrm{diag}\left\{ k_{xm}\right\} $. Произведение $[[\varepsilon]] \boldsymbol{E}_{y}$ соответствует свертке, так что Тёплицева матрица $[[\varepsilon]]$ (матрица, элементы которой зависят от разности индексов) составлена из коэффициентов Фурье периодической функции диэлектрической проницаемости: $[[\varepsilon]]_{mn} = \varepsilon_{m-n}$.
Из уравнений получаем систему дифференциальных уравнений \begin{equation} \dfrac{d^{2}\boldsymbol{E}_{y}}{dz^{2}}+\left(\omega^{2}\mu_{0} [[\varepsilon]] -K^{2}\right)\boldsymbol{E}_{y}=0 \end{equation} Подставляя решение в форме $\boldsymbol{E}_{x}=\boldsymbol{e}_{x}\exp\left(i\beta z\right)$ приходим к матрично-векторной задаче на собственные значения \begin{equation} \left(\omega^{2}\mu_{0} [[\varepsilon]] -K^{2}\right)\boldsymbol{e}_{y} = \beta^{2}\boldsymbol{e}_{y} \end{equation} которая устойчиво решается стандартными методами, реализованными во всех математических пакетах. После ее решения магнитное поле найденных мод рассчитывается из соотношения $\boldsymbol{h}_{xm}=-(\beta_m/\omega\mu_{0})\boldsymbol{e}_{ym}$. Полученное дисперсинное уравнение заменяет трансцендентное уравнение для модального метода.
В случае ТМ поляризации заметим, что компонента электрической индукции в области кристалла $D_x=\varepsilon(x)E_x$ является произведением двух функций, имеющих разрыв первого рода в одних и тех же точках, хотя сама по себе является непрерывной в этих точках. Производная же $D_x$ в точках разрыва не существует, поэтому сходимость конечных сумм в этих точках очень медленная, что будет опрелелять численную сходимость итогового метода. Чтобы улучшить сходимость метода применяют так называемые правила Ли, переписывая $D_x/\varepsilon(x)=E_x$, раскладывая это соотношение в ряд Фурье, обрывая ряды на конечных членах и записывая матрично-векторное соотношение на Фурье-амплитуды $\boldsymbol{D}_{x} = [[1/\varepsilon]]^{-1}\boldsymbol{E}_{x}$. Применяя этот подход, уравнения Максвелла для ТМ поляризации запишутся как \begin{align} \dfrac{d\boldsymbol{H}_{y}}{dz} &= i\omega [[1/\varepsilon]]^{-1}\boldsymbol{E}_{x} \\ K\boldsymbol{H}_{y} &= -\omega[[\varepsilon]] \boldsymbol{E}_{z} \\ \dfrac{d\boldsymbol{E}_{x}}{dz}-iK\boldsymbol{E}_{z} &= i\omega\mu_{0}\boldsymbol{H}_{y} \end{align} Далее поступая аналогично случаю ТЕ поляризации, можно получить следующее дисперсионное уравнение \begin{equation} \beta^{2} \boldsymbol{h}_{y} = \omega^{2}\mu_{0} [[1/\varepsilon]]^{-1} \left(\mathbb{I}-\dfrac{1}{\omega^{2}\mu_{0}}K [[\varepsilon]]^{-1} K\right) \boldsymbol{h}_{y} \end{equation} При этом $\boldsymbol{e}_{xm} = (\beta_m/\omega) [[1/\varepsilon]] \boldsymbol{h}_{ym}$.
Т-матрица плоской границы между кристаллом и однородной изотропной средой получается точно также, как и выше с тем отличием, что поле в кристалле теперь раскладывается по плоским волнам. Обозначая коэффициенты Фурье разлолжения собственного поля $q$-й моды как $f_{qm}$, $g_{qm}$ (это собственные вектора, найденные на предыдущем этапе), запишем разложение в явном виде: \begin{align} F\left(x,-0\right)=\sum_{m}e^{ik_{xm}x} \sum_{q}f_{yqm}\left(c_{q}^{e+}+c_{q}^{e-}\right) \\ G\left(x,-0\right)=\sum_{m}e^{ik_{xm}x} \sum_{q}g_{xqm}\left(c_{q}^{e+}-c_{q}^{e-}\right) \end{align} Условия сопряжения приводят к следующему виду Т-матрицы: \begin{equation} \left(\begin{array}{c} a_{m}^{+} \\ a_{m}^{-} \end{array}\right) = \dfrac{1}{2} \sum_{q} \left( \begin{array}{cc} f_{qm}-\dfrac{\omega\eta}{k_{zm}}g_{qm} & f_{qm}+\dfrac{\omega\eta}{k_{zm}}g_{qm} \\ f_{qm}+\dfrac{\omega\eta}{k_{zm}}g_{qm} & f_{qm}-\dfrac{\omega\eta}{k_{zm}}g_{qm} \end{array} \right) \left( \begin{array}{c} c_{q}^{+} \\ c_{q}^{-} \end{array} \right) \end{equation} Далее эта Т-матрица может быть использована для получения матрицы рассеяния слоя кристалла точно также, как это было сделано для модального метода.