{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5a3ba2b9",
   "metadata": {},
   "source": [
    "Лицензия MIT\n",
    "\n",
    "© Алексей Александрович Щербаков, 2024"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00ffec2d",
   "metadata": {},
   "source": [
    "# Лекция 3.3. Дифракционные решетки. Резонансные свойства"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab869914",
   "metadata": {},
   "source": [
    "## Резонансы и полюса матрицы рассеяния\n",
    "\n",
    "Во вводной части было показано, что моды фотонных структур можно изучать не только путем решения задачи на собственные значения для волнового уравнения, но и путем расчета полюсов матрицы рассеяния. Матрица рассеяния дает единый способ описания слоя, позволяя получить и дискретный и непрерывный спектр. Пусть на некоторую оптическую структуру падают волны, задаваемые вектором амплитуд $\\boldsymbol{a}_{inc}$, а рассеянное поле задается вектором амплитуд $\\boldsymbol{a}_{sca}$. Запишем связь между ними через обратную матрицу рассеяния как $S^{-1}\\boldsymbol{a}_{sca} = \\boldsymbol{a}_{inc}$. При нулевой правой части решениями этого уравнения будут являться собственные решения $\\boldsymbol{a}_{eig}$ уравнения Гельмгольца без источников, то есть, собственные числа будут определяться полюсами матрицы рассеяния, а собственные вектора - вычетами в этих полюсах.\n",
    "\n",
    "Полюса могут быть найдены одним из численных методов, рассмотренных в общем блоке. Если изучать резонансы в дифракционных решетках, можно использвать промежутоные результаты Фурье-модального метода, или же приближенно описать резонансы через интерференцию Блоховских мод. Рассмотрим ниже эти подходы."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e4ad999",
   "metadata": {},
   "source": [
    "## Моды решетки\n",
    "\n",
    "###  Расчет мод решетки в рамках ФММ\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",
    "\\begin{equation}\n",
    "    \\boldsymbol{b}^{+} = T_{21}^{(1)} \\mathcal{E}\\boldsymbol{c}^{-} + T_{22}^{(1)} \\boldsymbol{c}^{+}\n",
    "    \\Rightarrow \\boldsymbol{c}^{+} = \\left(T_{22}^{(1)}\\right)^{-1} \\boldsymbol{b}^{+} - \\left(T_{22}^{(1)}\\right)^{-1} T_{21}^{(1)} \\mathcal{E}\\boldsymbol{c}^{-} \\Rightarrow R_c^{(1)} = - \\left(T_{22}^{(1)}\\right)^{-1} T_{21}^{(1)}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "    \\boldsymbol{a}^{-} = T_{21}^{(2)} \\mathcal{E}\\boldsymbol{c}^{+} + T_{22}^{(2)} \\boldsymbol{c}^{-} \n",
    "    \\Rightarrow \\boldsymbol{c}^{-} = \\left(T_{22}^{(2)}\\right)^{-1} \\boldsymbol{a}^{-} - \\left(T_{22}^{(2)}\\right)^{-1} T_{21}^{(2)} \\mathcal{E}\\boldsymbol{c}^{+} \\Rightarrow R_c^{(2)} = - \\left(T_{22}^{(2)}\\right)^{-1} T_{21}^{(2)}\n",
    "\\end{equation}\n",
    "При осутствии падающего поля получаем условие на существование собственного решения:\n",
    "\\begin{equation}\n",
    "    \\left\\{ \\begin{array}{c} \\boldsymbol{c}^{+} = R_c^{(1)} \\mathcal{E}\\boldsymbol{c}^{-} \\\\ \\boldsymbol{c}^{-} = R_c^{(2)} \\mathcal{E}\\boldsymbol{c}^{+} \\end{array} \\right. \\Rightarrow \\det \\left( \\mathbb{I} - R_c^{(1)} \\mathcal{E} R_c^{(2)} \\mathcal{E} \\right) = 0\n",
    "\\end{equation}\n",
    "Решение этого уравнения даёт частоты собственных мод решетки.\n",
    "\n",
    "### Резонансы как интерференция блоховских мод\n",
    "\n",
    "Модальное представление также позволяет получить простое приближение и физически проинтерпретировать высокодобротные резонансы в диэлектрических решетках. Рассмотрим симметричную структуру, для которой $R_c^{(1)} = R_c^{(2)}$. Пусть период и контраст решетки таковы, что в ней существует небольшое число распространяющихся блоховских мод, для которых $\\beta_m^2 > 0$, $1\\leq m\\leq N_m$. Пренебрежем всеми затухающими модами и рассмотрим более подробно отражение мод на границах решетки.\n",
    "\n",
    "В одномодовом приближении, когда $N_m=1$, равенство нулю детерминанта приводит к условию\n",
    "\\begin{equation}\n",
    "    1-r_{11}^{2}\\exp\\left(2i\\beta_{1}h\\right)=0\\Rightarrow\\beta_{1}h+\\arg r_{11}=\\pi l\n",
    "\\end{equation}\n",
    "что показывает, что в данном приближении решетка ведет себя точно также, как и плоский однородный диэлектрический слой. Амплитуда уходящей нулевой гармоники\n",
    "\\begin{equation}\n",
    "    a_1^+ = t_{1} c_1^+ \\exp(i\\beta_1 h)\n",
    "\\end{equation}\n",
    "и $|t_{1}|\\neq 0$ для любого блоховского вектора.\n",
    "\n",
    "В двумодовом приближении, когда $N_m=2$ возникают ситуации, не наблюдающиеся для однородного слоя. Во-первых, в Г-точке первая мода является симметричной, а вторая - антисимметричной. Это значит, что интеграл перекрытия внешней плоской волны, распространяющейся нормально к решетке, равен нулю, а, следовательно, $r_{12}=r_{21}=0$ а также $r_{22}=1$ и $t_2 = 0$. Состояние с возбужденной второй модой, обладающее бесконечной добротностью в Г-точке, называется связанным состоянием в континууме, защищенным симметрией. Во-вторых, для $k_B\\neq0$ двумодовое приближение даёт четыре уравнения на амплитуды мод, которые можно переписать в той же форме, что и уравнения для резонатора Фабри-Перо:\n",
    "\\begin{equation}\n",
    "    1-\\left(r^{(12)}\\right)^{2}\\exp\\left(2i\\beta_{2}h\\right)=0\n",
    "\\end{equation}\n",
    "где\n",
    "\\begin{equation}\n",
    "    r^{(12)}=\\dfrac{r_{22}+\\alpha r_{11}r_{21}r_{12}\\exp\\left(2i\\beta_{1}h\\right)}{1-\\alpha r_{21}r_{12}\\exp\\left[i\\left(\\beta_{1}+\\beta_{2}\\right)h\\right]}, \\; \\;  \\alpha=\\dfrac{1}{1-r_{11}^{2}\\exp\\left(2i\\beta_{1}h\\right)}\n",
    "\\end{equation}\n",
    "Аналогично, амплитуду уходящей плоской волны можно выразить через эффективный коэффициент пропускания:\n",
    "\\begin{equation}\n",
    "    b^{t} = t_{1}a_{1}^{+} \\exp\\left(i\\beta_{1}h\\right) + t_{2}a_{2}^{+}\\exp\\left(i\\beta_{2}h\\right) = t^{(12)}a_{2}^{+}\\exp\\left(i\\beta_{2}h\\right)\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "    t^{(12)}=t_{2}+\\alpha t_{1}r_{21}\\exp\\left(i\\beta_{1}h\\right)\\left[r^{(12)}\\exp\\left(i\\beta_{2}h\\right)+r_{11}\\exp\\left(i\\beta_{1}h\\right)\\right]\n",
    "\\end{equation}\n",
    "который уже может полностью обнулится для некоторого значения блоховского вектора. Данную ситуацию можно интерпретировать как деструктивную интерференцию полей мод в пространстве, окружающем решетку. В точке, где обнуляется эффекивный коэффициент пропускания, добротность этого состояния также равна бесконечности, а состояние называется случайныи или интерференционным связанным состоянием в континууме, или ССК Фридриха-Винтгена."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15c9db4f",
   "metadata": {
    "vscode": {
     "languageId": "plaintext"
    }
   },
   "source": [
    "## Связь полюсного разложения и метода связанных волн\n",
    "\n",
    "Рассмотрим модели описания резонансного поведения волноводных решеток при дифракции монохроматической плоской волны, падающей на решетку под разными углами. В рамках подхода, который описывает резонансные свойства волноводных решеток через полюса мероморфных функций проекции волнового вектора в плоскости решетки, основными параметрами для построения кривой коэффициента отражения являтеся внерезонансный коэффициент отражения, который в некотором диапазоне изменения аргумента резонансной функции вблизи резонанса можно считать постоянным, координаты полюса на комплесной плоскости и амплитуда полюса. С другой стороны существует более физически интуитивный подход описания резонансных свойств волноводных решеток через формализм связанных мод. При этом каждая мода описывается постоянной распространения $\\beta$, коэффициентом излучения $\\alpha$, который характеризует потери энергии моды за счет наличия решетки (в плоском волноводе он равен нулю), и коэффициентом связи этой моды с падающей на решетку плоской волной $\\kappa$. Условия возбуждения (сохранение импульса в проекции на плоскость решетки) описываются с помощью оптогеометрического параметра\n",
    "\\begin{equation}\n",
    "    k_x  = k_x^{inc} \\pm K_g = k_0 n_c \\sin\\theta_{inc} \\pm K_g\n",
    "\\end{equation}\n",
    "где $k_0=2\\pi/\\lambda$, и знак $\\pm$ соответствует связи падающего излучения с модой, которая распространяется в том же направлении, что и падающая волна, либо в противоположном.\n",
    "\n",
    "Ценность описания резонансного поведения с помощью формализма связанных мод заключается в том, что он не только дает физически интуитивное описания возникновения резонанса, но и позволяет решать обратную задачу проектирования резонансных устройств на интуитивном уровне. Проблемой этого подхода является то, что параметры $\\alpha$, $\\beta$ и $\\kappa$ невозможно измерить напрямую. Рассмотрим, как можно получить эти параметры анализируя зависимость коэффициента отражения от проекции волнового вектора $r(k_x)$.\n",
    "\n",
    "<figure>\n",
    "    <center>\n",
    "    <img src=\"../pic/3-3_coupled_mode.png\" width=\"600\" align=\"center\"/>\n",
    "        <figcaption>Схема возбуждения решеточной моды, описывающая резонансное отражение</figcaption>\n",
    "    </center>\n",
    "</figure>\n",
    "\n",
    "Необходимое условие возбуждение моды при коллинеарном падении плоской волны с помощью \"$-1$\"-го дифракционного порядка в k-пространстве записывается как\n",
    "\\begin{equation}\n",
    "    \\beta = -n_m k_0 = k_x^{inc} - K_g\n",
    "\\end{equation}\n",
    "Обозначим показатели преломления подложки, волновода и покрытия как $n_s$, $n_g$ и $n_c$, так что $n_g>n_s>n_c$, $n_g\\equiv n_m$. В отражение вносят вклад следующие компоненты. Во-первых, это френелевское отражение, которое обозначим символом $r_0$. В-вторых, это вклад $r_g$, обусловленный связанной модой, которая переизлучается за счет присутствия решетки. Запишем зависиомость $r(k_x)$ с помощью вклада полюсного слагаемого\n",
    "\\begin{equation}\n",
    "    r(k_x) = r_0 + \\frac{a_p}{k_x - k_p}\n",
    "\\end{equation}\n",
    "Переменная $k_x$ зависит от параметров задачи $\\lambda$ и $\\Lambda$, а также от угла падения $\\theta_{inc}$, так что экспериментально и с помощью моделирования можно напряму получать $r(k_x)$ как функцию действительного переменного путем изменения угла падения.\n",
    "\n",
    "Амплитуда волноводной моды $C_g(x)$ удовлетворяет дифференциальному уравнению, означающему то, что скорость изменения амплитуды моды обусловлена возбуждением внешним полем и излучением за счет решетки:\n",
    "\\begin{equation}\n",
    "    \\frac{\\partial C}{\\partial x} = \\kappa q(x) \\exp[i(k_x-\\beta)x] - \\alpha C_g\n",
    "\\end{equation}\n",
    "где $q(x)=1$ - профиль поля бескончной плоской волны, падающей на бесконечную решетку. Разность $k_x-\\beta$ показывает расстройку от идеального синхронизма возбуждения моды. Стационарное решение этого уравнения в пределе условия синхронизма\n",
    "\\begin{equation}\n",
    "    C\\left(x\\right)=c\\left(x\\right)e^{-\\alpha x}\\Rightarrow\\dfrac{dc}{dx}=\\kappa e^{\\alpha x+i\\left(k_x-\\beta\\right)x}\\Rightarrow C\\left(x\\right)=\\frac{-i\\kappa}{k_x-\\beta-i\\alpha}\n",
    "\\end{equation}\n",
    "Полученная амплитуда моды должна быть линейно связана с решеточным вкладом в коэффициент отражения\n",
    "\\begin{equation}\n",
    "    r_g(k_x) = \\eta C = \\frac{-i\\eta\\kappa}{k_x-(\\beta+i\\alpha)}\n",
    "\\end{equation}\n",
    "Физический смысл коэффиента $\\eta$ состоит в том, что он определяет дифракционные потери моды в покрытие. Сравнивая полученное выражение с полюсным разложением, имеем\n",
    "\\begin{split}\n",
    "    a_p = -i\\eta\\kappa \\\\\n",
    "    k_p = \\beta+i\\alpha\n",
    "\\end{split}\n",
    "\n",
    "Рассмотрим графическую интерпретацию полученных соотношений на комплексной плоскости коэффициента отражения. При изменении параметра $k_x$ от  $-\\infty$ до $+\\infty$ коэффициент отражения описывает окружность на комплесной плоскости. Это следует из того, что легко показать, что\n",
    "\\begin{equation}\n",
    "    \\left( \\Re e[r_g] + \\frac{\\Im m[a_p]}{2\\alpha} \\right)^2 + \\left( \\Im e[r_g] - \\frac{\\Re m[a_p]}{2\\alpha} \\right)^2 = \\left( \\frac{|a_p|}{2\\alpha} \\right)^2\n",
    "\\end{equation}\n",
    "Центр этой окружности - точка $z_0 = ia_p/2\\alpha$, а радиус равен $\\rho = a_p|/2\\alpha$. Положение точки на окружности определяется аргументом\n",
    "\\begin{equation}\n",
    "    \\mathrm{Arg} [r_g(k_x) - z_0] = \\varphi_p - \\pi/2 - 2\\arctan\\left(\\frac{k_x-\\beta}{\\alpha}\\right)\n",
    "\\end{equation}\n",
    "Вне резонанса коэффициент отржения $r_g$ становится равным нулю, поэтому окружность пересекает 0.\n",
    "\n",
    "<figure>\n",
    "    <center>\n",
    "    <img src=\"../pic/3-3_complex_reflection.png\" width=\"500\" align=\"center\"/>\n",
    "        <figcaption>Комплексный коэффициент отражения с учетом решеточного вклада</figcaption>\n",
    "    </center>\n",
    "</figure>\n",
    "\n",
    "Полный коэффициенто тражения складывается из отражения Фабри-Перо и отражения, обусловленного волноводной решеточной модой, поэтому окружность полного коэффициента отражения смещается на $r_0$ в точку\n",
    "\\begin{equation}\n",
    "    z_c = r_0 + ia_p/2\\alpha\n",
    "\\end{equation}\n",
    "Если параметры решетки подобраны так, что модуль коэффиента отражения в резонансе равен 1, то окружность $r(k_x)$ касается единичной окружности в точке $r_{max}$. При этом, очевидно, точки максимального $r_{max}$ и минимального $r_{min}$ коэффицента отражения, а также центр $z_c$ лежат на прямой, проходящей через ноль комплексной плоскости. Условие резонансного возбуждения моды $k_x=\\beta$ дает выражение для коэффициента отражения в этой точке\n",
    "\\begin{equation*}\n",
    "    r(k_x=\\beta) = r_{beta} = r_0 + \\frac{ia_p}{\\alpha}\n",
    "\\end{equation*}\n",
    "Отсюда видно, что $r_{\\beta}$, $r_0$ и $z_c$ лежать на одной прямой. Тогда\n",
    "\\begin{equation*}\n",
    "    z_c = \\frac{1}{2}(r_0 + r_{\\beta}) = \\frac{1}{2}(r_{max} + r_{min}) \\Rightarrow r_{\\beta} = r_{max} + r_{min} - r_0\n",
    "\\end{equation*}\n",
    "то есть, $r_{\\beta}$ оказывается выражен через измеримые или легко рассчитываемые величины.\n",
    "\n",
    "Используя полученные соотношения, можно показать, что модуль коэффицента отражения\n",
    "\\begin{equation*}\n",
    "    |r(k_x)|^2 = \\frac{|r_0|^2(k_x-\\beta)^2 + |r_\\beta|^2\\alpha^2 + 2\\alpha(k_x-\\beta)\\zeta}{(k_x-\\beta)^2 + \\alpha^2}\n",
    "\\end{equation*}\n",
    "где $\\zeta = \\sqrt{(|r_{max}|^2-|r_0|^2)/(|r_0|^2-|r_{min}|^2)}$. Из условия $d(|r(k_x)|^2)/dk_x = 0$ можно найти проекции волнового вектора, соответствующие максимуму и минимуму отражения, окуда получаются явные выражения для коэффиента дифракционных потерь $\\alpha$ и константы распространения моды $\\beta$:\n",
    "\\begin{equation*}\n",
    "    \\begin{array}{l} k_{x,max}-\\beta = \\alpha/\\zeta \\\\ k_{x,min}-\\beta = \\alpha\\zeta \\end{array} \\Rightarrow \\alpha = \\zeta \\frac{k_{x,max}-k_{x,min}}{1+\\zeta^2}\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "    \\beta = \\frac{k_{x,max}+k_{x,min}}{2} + \\frac{k_{x,max}-k_{x,min}}{2} \\frac{\\zeta^2-1}{\\zeta^2+1}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5ba84dc",
   "metadata": {},
   "source": [
    "#### Литература\n",
    "\n",
    "1. 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",
    "2. P. Lalanne, J. P. Hugonin, and P. Chavel, [Optical Properties of Deep Lamellar Gratings: A Coupled Bloch-Mode Insight](https://opg.optica.org/jlt/abstract.cfm?uri=JLT-24-6-2442), J. Lightwave Technol. 24, 2442- (2006)\n",
    "3. A. Ovcharenko, C. Blanchard, J. P. Hugonin, and C. Sauvar, [Bound states in the continuum in symmetric and asymmetric photonic crystal slabs](https://doi.org/10.1103/PhysRevB.101.155303), Phys. Rev. B 101, 155303 (2020)\n",
    "4. G. Quaranta, G. Basset, O. J. F. Martin, B. Gallinet, [Recent Advances in Resonant Waveguide Gratings](https://doi.org/10.1002/lpor.201800017) Laser & Photonics Reviews 12, 1800017 (2018)\n",
    "5. D. Pietroy, A. V. Tishchenko, M. Flury, and O. Parriaux, [Bridging pole and coupled wave formalisms for grating waveguide resonance analysis and design synthesis](https://opg.optica.org/oe/fulltext.cfm?uri=oe-15-15-9831&id=140070), Opt. Express 15, 9831-9842 (2007)\n",
    "6. D. Rosenblatt, A. Sharon and A. A. Friesem, [Resonant grating waveguide structures](https://doi.org/10.1109/3.641320), in IEEE Journal of Quantum Electronics, vol. 33, no. 11, pp. 2038-2059, Nov. 1997"
   ]
  }
 ],
 "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
}