{ "cells": [ { "cell_type": "markdown", "id": "5a3ba2b9", "metadata": {}, "source": [ "Лицензия MIT\n", "\n", "© Алексей Александрович Щербаков, 2024" ] }, { "cell_type": "markdown", "id": "00ffec2d", "metadata": {}, "source": [ "# Лекция 3.2. Дифракционные решетки. Модальные методы." ] }, { "cell_type": "markdown", "id": "ab869914", "metadata": {}, "source": [ "Модальные методы можно применять для решения различных задач с различной геометрией. Основная идея заключается в том, чтобы рассчитать полный модальный базис в различных областях пространства, а затем произвести сшивку базисов на границах соседних областей при помощи условий сопряжения. При этом, ввиду взаимной неортогональности базисов в разных подпространствах, возникают интегралы перекрытия различных мод. Здесь эта идея будет продемонстрирована на примере слоя одномерного фотонного кристалла." ] }, { "cell_type": "markdown", "id": "219c6483", "metadata": {}, "source": [ "## T-матрица границы фотонного кристалла\n", "\n", "Воспользуемся результатами модального разложения для одномерного бесконечного фотонного кристалла. Пусть теперь ось $X$ декартовых координат задает направление периодичности. Рассмотрим плоскую границу $z=0$ такую, что полупространство ниже границы при $z<0$ заполнено фотонным кристаллом, а полупространство $z>0$ заполнено однородной изотропной средой с параметрами $\\varepsilon$, $\\mu$. Зафиксируем блоховский волновой вектор $k_B \\equiv k_{x0}$. Разложение поля в верхнем полупространстве по плоским волнам имеет вид\n", "\\begin{align*}\n", " 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} \\\\\n", " 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}\n", "\\end{align*}\n", "где $k_{xm} = k_{x0} + 2\\pi m/\\Lambda$, $k_{zm} = \\sqrt{\\omega^2\\varepsilon\\mu - k_{xm}^2}$. Счетный набор плоских волн в разложении возникает вследствие заданной фотонным кристаллом периодичности всего пространства. В области $z<0$ воспользуемся полученным модальным разложением:\n", "\\begin{align*}\n", " 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) \\\\\n", " 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)\n", "\\end{align*}\n", "Непрерывность тангенциальных компонент полей на границе $z=0$ даёт:\n", "\\begin{align*}\n", " \\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} \\\\\n", " \\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}\n", "\\end{align*}\n", "Эту систему можно разрешить, например, для коэффициентов разложения поля в однородной среде. Воспользуемся ортогональностью плоских волн - умножая обе части уравнений на $\\exp(-ik_{xp}x)$ и интегрируя по периоду, получаем\n", "\\begin{align}\n", " 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) \\\\\n", " \\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)\n", "\\end{align}\n", "где введены обозначения интегралов перекрытия $\\mathcal{I}_{nm}^{(1,2)}$. Запишем связь между наборами коэффициентов разложения полей по обе стороны границы раздела полупространств в форме Т-матрицы:\n", "\\begin{equation}\n", " \\left(\\begin{array}{c}\n", "a_{n}^{+}\\\\\n", "a_{n}^{-}\n", "\\end{array}\\right)=\\sum_{m}\\left(\\begin{array}{cc}\n", "T_{11,nm} & T_{12,nm}\\\\\n", "T_{21,nm} & T_{22,nm}\n", "\\end{array}\\right)\\left(\\begin{array}{c}\n", "c_{n}^{+}\\\\\n", "c_{n}^{-}\n", "\\end{array}\\right) = \\dfrac{1}{2}\\left(\\begin{array}{cc}\n", "\\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)}\\\\\n", "\\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)}\n", "\\end{array}\\right) \\left(\\begin{array}{c}\n", "c_{n}^{+}\\\\\n", "c_{n}^{-}\n", "\\end{array}\\right)\n", "\\end{equation}\n", "В случае простой ячейки фотонного кристалла, состоящей из двух однородных слоёв, интегралы перекрытия легко находятся аналитически:\n", "\\begin{align}\n", " \\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) \\\\\n", " &+ \\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)\n", "\\end{align}\n", "\n", "\\begin{align}\n", " \\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) \\\\\n", " &+ \\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)\n", "\\end{align}" ] }, { "cell_type": "markdown", "id": "65c62167", "metadata": {}, "source": [ "## S-матрица слоя фотонного кристалла\n", "\n", "Далее перейдем к случаю, когда фотонный кристалл занимает область пространства в пределах плоского слоя толщиной $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}$ - константы распространения мод, являющиеся решениями дсперсионного уравнения в бесконечном фотонном кристалле.\n", "\n", "Воспользуемся выведенными Т-матрицами границ радела и запишем их для верхней и нижней границ рассматриваемого слоя:\n", "\\begin{equation}\n", " \\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)\n", "\\end{equation}\n", "\n", "\\begin{equation}\n", " \\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)\n", "\\end{equation}\n", "\n", "Эти соотношения можно переписать в форме S-матрицы:\n", "\\begin{equation}\n", " \\left(\\begin{array}{c} \\boldsymbol{b}^{-}\\\\ \\boldsymbol{a}^{+} \\end{array}\\right) = S \\left(\\begin{array}{c} \\boldsymbol{b}^{+}\\\\ \\boldsymbol{a}^{-} \\end{array}\\right)\n", "\\end{equation}\n", "\n", "где в явном виде\n", "\\begin{equation}\n", " 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)}\\\\\n", "T_{22}^{(2)} & T_{21}^{(2)}\\mathcal{E} \\end{array}\\right)^{-1}\n", "\\end{equation}\n", "Использование более простой для аналитических вычислений T-матрицы в данном случае приводит к вычислительной неустойчивости и расходимости метода, посколько в матрице переноса появляются множители $\\exp\\left(-i\\beta_{m}|\\Delta z|\\right)$, которые могут оказаться очень большими для сильно затухающих волн, что приводит к быстрому накоплению ошибки численных расчетов. Затухающие волны же необходимо учитывать при использовании данного численного метода расчета дифракции." ] }, { "cell_type": "markdown", "id": "f66b57a6", "metadata": {}, "source": [ "## Фурье-модальный метод\n", "\n", "Основными недостатками модального метода, описанного выше, является, во-первых, необходимость решения трансцендентного дисперсионного уравнения, при этом не пропуская его последовательные корни, а, во-вторых, сложность масштабирования подхода на двумерные фотонные кристаллы. Оба недостатка устраняются в рамках так Фурье-модального метода. Последовательность шагов остаётся той же, что при формулировке модального метода, а именно, решение задачи на собственные значения в бесконечном фотонном кристалле, нахождение Т-матрицы плоской границы, и нахождение матрицы рассеяния плоского слоя кристалла. Третий шаг идентичен уже описанному выше, поэтому здесь приведем детали первых двух шагов.\n", "\n", "### Дисперсионное уравнение\n", "\n", "Для всех компонент полей кроме теоремы Блоха запишем разложение в ряд Фурье периодической части волновой функции:\n", "\\begin{equation}\n", " 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)\n", "\\end{equation}\n", "Подставляя подобные разложения в уравнения Максвелла для двумерной задачи и ТЕ-поляризованного поля, получим\n", "\\begin{align*}\n", " -\\frac{d\\boldsymbol{E}_{y}}{dz} &= i\\omega\\mu_{0}\\boldsymbol{H}_{x} \\\\\n", " K\\boldsymbol{E}_{y} &= \\omega\\mu_{0}\\boldsymbol{H}_{z} \\\\\n", " \\dfrac{d\\boldsymbol{H}_{x}}{dz} - iK\\boldsymbol{H}_{z} &= -i\\omega [[\\varepsilon]] \\boldsymbol{E}_{y}\n", "\\end{align*}\n", "где векторы $\\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}$.\n", "\n", "Из уравнений получаем систему дифференциальных уравнений\n", "\\begin{equation}\n", " \\dfrac{d^{2}\\boldsymbol{E}_{y}}{dz^{2}}+\\left(\\omega^{2}\\mu_{0} [[\\varepsilon]] -K^{2}\\right)\\boldsymbol{E}_{y}=0\n", "\\end{equation}\n", "Подставляя решение в форме $\\boldsymbol{E}_{x}=\\boldsymbol{e}_{x}\\exp\\left(i\\beta z\\right)$ приходим к матрично-векторной задаче на собственные значения\n", "\\begin{equation}\n", " \\left(\\omega^{2}\\mu_{0} [[\\varepsilon]] -K^{2}\\right)\\boldsymbol{e}_{y} = \\beta^{2}\\boldsymbol{e}_{y}\n", "\\end{equation}\n", "которая устойчиво решается стандартными методами, реализованными во всех математических пакетах. После ее решения магнитное поле найденных мод рассчитывается из соотношения $\\boldsymbol{h}_{xm}=-(\\beta_m/\\omega\\mu_{0})\\boldsymbol{e}_{ym}$. Полученное дисперсинное уравнение заменяет трансцендентное уравнение для модального метода.\n", "\n", "В случае ТМ поляризации заметим, что компонента электрической индукции в области кристалла $D_x=\\varepsilon(x)E_x$ является произведением двух функций, имеющих разрыв первого рода в одних и тех же точках, хотя сама по себе является непрерывной в этих точках. Производная же $D_x$ в точках разрыва не существует, поэтому сходимость конечных сумм в этих точках очень медленная, что будет опрелелять численную сходимость итогового метода. Чтобы улучшить сходимость метода применяют так называемые правила Ли, переписывая $D_x/\\varepsilon(x)=E_x$, раскладывая это соотношение в ряд Фурье, обрывая ряды на конечных членах и записывая матрично-векторное соотношение на Фурье-амплитуды $\\boldsymbol{D}_{x} = [[1/\\varepsilon]]^{-1}\\boldsymbol{E}_{x}$. Применяя этот подход, уравнения Максвелла для ТМ поляризации запишутся как\n", "\\begin{align}\n", " \\dfrac{d\\boldsymbol{H}_{y}}{dz} &= i\\omega [[1/\\varepsilon]]^{-1}\\boldsymbol{E}_{x} \\\\\n", " K\\boldsymbol{H}_{y} &= -\\omega[[\\varepsilon]] \\boldsymbol{E}_{z} \\\\\n", " \\dfrac{d\\boldsymbol{E}_{x}}{dz}-iK\\boldsymbol{E}_{z} &= i\\omega\\mu_{0}\\boldsymbol{H}_{y}\n", "\\end{align}\n", "Далее поступая аналогично случаю ТЕ поляризации, можно получить следующее дисперсионное уравнение\n", "\\begin{equation}\n", " \\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}\n", "\\end{equation}\n", "При этом $\\boldsymbol{e}_{xm} = (\\beta_m/\\omega) [[1/\\varepsilon]] \\boldsymbol{h}_{ym}$.\n", "\n", "### T-матрица плоской границы\n", "\n", "Т-матрица плоской границы между кристаллом и однородной изотропной средой получается точно также, как и выше с тем отличием, что поле в кристалле теперь раскладывается по плоским волнам. Обозначая коэффициенты Фурье разлолжения собственного поля $q$-й моды как $f_{qm}$, $g_{qm}$ (это собственные вектора, найденные на предыдущем этапе), запишем разложение в явном виде:\n", "\\begin{align}\n", " F\\left(x,-0\\right)=\\sum_{m}e^{ik_{xm}x} \\sum_{q}f_{yqm}\\left(c_{q}^{e+}+c_{q}^{e-}\\right) \\\\\n", " G\\left(x,-0\\right)=\\sum_{m}e^{ik_{xm}x} \\sum_{q}g_{xqm}\\left(c_{q}^{e+}-c_{q}^{e-}\\right)\n", "\\end{align}\n", "Условия сопряжения приводят к следующему виду Т-матрицы:\n", "\\begin{equation}\n", " \\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} \\\\\n", "f_{qm}+\\dfrac{\\omega\\eta}{k_{zm}}g_{qm} & f_{qm}-\\dfrac{\\omega\\eta}{k_{zm}}g_{qm}\n", "\\end{array} \\right) \\left( \\begin{array}{c} c_{q}^{+} \\\\ c_{q}^{-} \\end{array} \\right)\n", "\\end{equation}\n", "Далее эта Т-матрица может быть использована для получения матрицы рассеяния слоя кристалла точно также, как это было сделано для модального метода." ] }, { "cell_type": "markdown", "id": "fdc22726", "metadata": {}, "source": [ "#### Литература\n", "\n", "1. I. C. Botten, M. S. Craig, R. C. McPhedran, J. L. Adams, and J. R. Andrewartha, [The Dielectric Lamellar Diffraction Grating](https://doi.org/10.1080/713820571), Optica Acta: International Journal of Optics, 28(3), 413–428 (1981)\n", "2. A. V. Tishchenko, [Phenomenological representation of deep and high contrast lamellar gratings by means of the modal method](https://doi.org/10.1007/s11082-005-1188-2), Opt. Quant. Electron. 37, 309–330 (2005)\n", "3. [Gratings: Theory and Numeric Applications, Second Revisited Edition](https://www.fresnel.fr/files/gratings/Second-Edition/), ed. by E. Popov, Institut Fresnel, AMU, CNRS, ECM (2014)\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" } }, "nbformat": 4, "nbformat_minor": 5 }