#!/usr/bin/env python
# coding: utf-8

# Лицензия MIT
# 
# © Алексей Александрович Щербаков, 2024

# # Лекция 3.3. Дифракционные решетки. Резонансные свойства

# ## Резонансы и полюса матрицы рассеяния
# 
# Во вводной части было показано, что моды фотонных структур можно изучать не только путем решения задачи на собственные значения для волнового уравнения, но и путем расчета полюсов матрицы рассеяния. Матрица рассеяния дает единый способ описания слоя, позволяя получить и дискретный и непрерывный спектр. Пусть на некоторую оптическую структуру падают волны, задаваемые вектором амплитуд $\boldsymbol{a}_{inc}$, а рассеянное поле задается вектором амплитуд $\boldsymbol{a}_{sca}$. Запишем связь между ними через обратную матрицу рассеяния как $S^{-1}\boldsymbol{a}_{sca} = \boldsymbol{a}_{inc}$. При нулевой правой части решениями этого уравнения будут являться собственные решения $\boldsymbol{a}_{eig}$ уравнения Гельмгольца без источников, то есть, собственные числа будут определяться полюсами матрицы рассеяния, а собственные вектора - вычетами в этих полюсах.
# 
# Полюса могут быть найдены одним из численных методов, рассмотренных в общем блоке. Если изучать резонансы в дифракционных решетках, можно использвать промежутоные результаты Фурье-модального метода, или же приближенно описать резонансы через интерференцию Блоховских мод. Рассмотрим ниже эти подходы.

# ## Моды решетки
# 
# ###  Расчет мод решетки в рамках ФММ
# 
# В рамках Фурье-модального метода в явном виде были получены Т-матрицы границ раздела слоя фотонного кристалла с однородными средами:
# \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}
# Из этих уравнений можно получить матрицы отражения для модальных коэффициентов на границах слоя:
# \begin{equation}
#     \boldsymbol{b}^{+} = T_{21}^{(1)} \mathcal{E}\boldsymbol{c}^{-} + T_{22}^{(1)} \boldsymbol{c}^{+}
#     \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)}
# \end{equation}
# 
# \begin{equation}
#     \boldsymbol{a}^{-} = T_{21}^{(2)} \mathcal{E}\boldsymbol{c}^{+} + T_{22}^{(2)} \boldsymbol{c}^{-} 
#     \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)}
# \end{equation}
# При осутствии падающего поля получаем условие на существование собственного решения:
# \begin{equation}
#     \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
# \end{equation}
# Решение этого уравнения даёт частоты собственных мод решетки.
# 
# ### Резонансы как интерференция блоховских мод
# 
# Модальное представление также позволяет получить простое приближение и физически проинтерпретировать высокодобротные резонансы в диэлектрических решетках. Рассмотрим симметричную структуру, для которой $R_c^{(1)} = R_c^{(2)}$. Пусть период и контраст решетки таковы, что в ней существует небольшое число распространяющихся блоховских мод, для которых $\beta_m^2 > 0$, $1\leq m\leq N_m$. Пренебрежем всеми затухающими модами и рассмотрим более подробно отражение мод на границах решетки.
# 
# В одномодовом приближении, когда $N_m=1$, равенство нулю детерминанта приводит к условию
# \begin{equation}
#     1-r_{11}^{2}\exp\left(2i\beta_{1}h\right)=0\Rightarrow\beta_{1}h+\arg r_{11}=\pi l
# \end{equation}
# что показывает, что в данном приближении решетка ведет себя точно также, как и плоский однородный диэлектрический слой. Амплитуда уходящей нулевой гармоники
# \begin{equation}
#     a_1^+ = t_{1} c_1^+ \exp(i\beta_1 h)
# \end{equation}
# и $|t_{1}|\neq 0$ для любого блоховского вектора.
# 
# В двумодовом приближении, когда $N_m=2$ возникают ситуации, не наблюдающиеся для однородного слоя. Во-первых, в Г-точке первая мода является симметричной, а вторая - антисимметричной. Это значит, что интеграл перекрытия внешней плоской волны, распространяющейся нормально к решетке, равен нулю, а, следовательно, $r_{12}=r_{21}=0$ а также $r_{22}=1$ и $t_2 = 0$. Состояние с возбужденной второй модой, обладающее бесконечной добротностью в Г-точке, называется связанным состоянием в континууме, защищенным симметрией. Во-вторых, для $k_B\neq0$ двумодовое приближение даёт четыре уравнения на амплитуды мод, которые можно переписать в той же форме, что и уравнения для резонатора Фабри-Перо:
# \begin{equation}
#     1-\left(r^{(12)}\right)^{2}\exp\left(2i\beta_{2}h\right)=0
# \end{equation}
# где
# \begin{equation}
#     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)}
# \end{equation}
# Аналогично, амплитуду уходящей плоской волны можно выразить через эффективный коэффициент пропускания:
# \begin{equation}
#     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)
# \end{equation}
# 
# \begin{equation}
#     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]
# \end{equation}
# который уже может полностью обнулится для некоторого значения блоховского вектора. Данную ситуацию можно интерпретировать как деструктивную интерференцию полей мод в пространстве, окружающем решетку. В точке, где обнуляется эффекивный коэффициент пропускания, добротность этого состояния также равна бесконечности, а состояние называется случайныи или интерференционным связанным состоянием в континууме, или ССК Фридриха-Винтгена.

# ## Связь полюсного разложения и метода связанных волн
# 
# Рассмотрим модели описания резонансного поведения волноводных решеток при дифракции монохроматической плоской волны, падающей на решетку под разными углами. В рамках подхода, который описывает резонансные свойства волноводных решеток через полюса мероморфных функций проекции волнового вектора в плоскости решетки, основными параметрами для построения кривой коэффициента отражения являтеся внерезонансный коэффициент отражения, который в некотором диапазоне изменения аргумента резонансной функции вблизи резонанса можно считать постоянным, координаты полюса на комплесной плоскости и амплитуда полюса. С другой стороны существует более физически интуитивный подход описания резонансных свойств волноводных решеток через формализм связанных мод. При этом каждая мода описывается постоянной распространения $\beta$, коэффициентом излучения $\alpha$, который характеризует потери энергии моды за счет наличия решетки (в плоском волноводе он равен нулю), и коэффициентом связи этой моды с падающей на решетку плоской волной $\kappa$. Условия возбуждения (сохранение импульса в проекции на плоскость решетки) описываются с помощью оптогеометрического параметра
# \begin{equation}
#     k_x  = k_x^{inc} \pm K_g = k_0 n_c \sin\theta_{inc} \pm K_g
# \end{equation}
# где $k_0=2\pi/\lambda$, и знак $\pm$ соответствует связи падающего излучения с модой, которая распространяется в том же направлении, что и падающая волна, либо в противоположном.
# 
# Ценность описания резонансного поведения с помощью формализма связанных мод заключается в том, что он не только дает физически интуитивное описания возникновения резонанса, но и позволяет решать обратную задачу проектирования резонансных устройств на интуитивном уровне. Проблемой этого подхода является то, что параметры $\alpha$, $\beta$ и $\kappa$ невозможно измерить напрямую. Рассмотрим, как можно получить эти параметры анализируя зависимость коэффициента отражения от проекции волнового вектора $r(k_x)$.
# 
# <figure>
#     <center>
#     <img src="../pic/3-3_coupled_mode.png" width="600" align="center"/>
#         <figcaption>Схема возбуждения решеточной моды, описывающая резонансное отражение</figcaption>
#     </center>
# </figure>
# 
# Необходимое условие возбуждение моды при коллинеарном падении плоской волны с помощью "$-1$"-го дифракционного порядка в k-пространстве записывается как
# \begin{equation}
#     \beta = -n_m k_0 = k_x^{inc} - K_g
# \end{equation}
# Обозначим показатели преломления подложки, волновода и покрытия как $n_s$, $n_g$ и $n_c$, так что $n_g>n_s>n_c$, $n_g\equiv n_m$. В отражение вносят вклад следующие компоненты. Во-первых, это френелевское отражение, которое обозначим символом $r_0$. В-вторых, это вклад $r_g$, обусловленный связанной модой, которая переизлучается за счет присутствия решетки. Запишем зависиомость $r(k_x)$ с помощью вклада полюсного слагаемого
# \begin{equation}
#     r(k_x) = r_0 + \frac{a_p}{k_x - k_p}
# \end{equation}
# Переменная $k_x$ зависит от параметров задачи $\lambda$ и $\Lambda$, а также от угла падения $\theta_{inc}$, так что экспериментально и с помощью моделирования можно напряму получать $r(k_x)$ как функцию действительного переменного путем изменения угла падения.
# 
# Амплитуда волноводной моды $C_g(x)$ удовлетворяет дифференциальному уравнению, означающему то, что скорость изменения амплитуды моды обусловлена возбуждением внешним полем и излучением за счет решетки:
# \begin{equation}
#     \frac{\partial C}{\partial x} = \kappa q(x) \exp[i(k_x-\beta)x] - \alpha C_g
# \end{equation}
# где $q(x)=1$ - профиль поля бескончной плоской волны, падающей на бесконечную решетку. Разность $k_x-\beta$ показывает расстройку от идеального синхронизма возбуждения моды. Стационарное решение этого уравнения в пределе условия синхронизма
# \begin{equation}
#     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}
# \end{equation}
# Полученная амплитуда моды должна быть линейно связана с решеточным вкладом в коэффициент отражения
# \begin{equation}
#     r_g(k_x) = \eta C = \frac{-i\eta\kappa}{k_x-(\beta+i\alpha)}
# \end{equation}
# Физический смысл коэффиента $\eta$ состоит в том, что он определяет дифракционные потери моды в покрытие. Сравнивая полученное выражение с полюсным разложением, имеем
# \begin{split}
#     a_p = -i\eta\kappa \\
#     k_p = \beta+i\alpha
# \end{split}
# 
# Рассмотрим графическую интерпретацию полученных соотношений на комплексной плоскости коэффициента отражения. При изменении параметра $k_x$ от  $-\infty$ до $+\infty$ коэффициент отражения описывает окружность на комплесной плоскости. Это следует из того, что легко показать, что
# \begin{equation}
#     \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
# \end{equation}
# Центр этой окружности - точка $z_0 = ia_p/2\alpha$, а радиус равен $\rho = a_p|/2\alpha$. Положение точки на окружности определяется аргументом
# \begin{equation}
#     \mathrm{Arg} [r_g(k_x) - z_0] = \varphi_p - \pi/2 - 2\arctan\left(\frac{k_x-\beta}{\alpha}\right)
# \end{equation}
# Вне резонанса коэффициент отржения $r_g$ становится равным нулю, поэтому окружность пересекает 0.
# 
# <figure>
#     <center>
#     <img src="../pic/3-3_complex_reflection.png" width="500" align="center"/>
#         <figcaption>Комплексный коэффициент отражения с учетом решеточного вклада</figcaption>
#     </center>
# </figure>
# 
# Полный коэффициенто тражения складывается из отражения Фабри-Перо и отражения, обусловленного волноводной решеточной модой, поэтому окружность полного коэффициента отражения смещается на $r_0$ в точку
# \begin{equation}
#     z_c = r_0 + ia_p/2\alpha
# \end{equation}
# Если параметры решетки подобраны так, что модуль коэффиента отражения в резонансе равен 1, то окружность $r(k_x)$ касается единичной окружности в точке $r_{max}$. При этом, очевидно, точки максимального $r_{max}$ и минимального $r_{min}$ коэффицента отражения, а также центр $z_c$ лежат на прямой, проходящей через ноль комплексной плоскости. Условие резонансного возбуждения моды $k_x=\beta$ дает выражение для коэффициента отражения в этой точке
# \begin{equation*}
#     r(k_x=\beta) = r_{beta} = r_0 + \frac{ia_p}{\alpha}
# \end{equation*}
# Отсюда видно, что $r_{\beta}$, $r_0$ и $z_c$ лежать на одной прямой. Тогда
# \begin{equation*}
#     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
# \end{equation*}
# то есть, $r_{\beta}$ оказывается выражен через измеримые или легко рассчитываемые величины.
# 
# Используя полученные соотношения, можно показать, что модуль коэффицента отражения
# \begin{equation*}
#     |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}
# \end{equation*}
# где $\zeta = \sqrt{(|r_{max}|^2-|r_0|^2)/(|r_0|^2-|r_{min}|^2)}$. Из условия $d(|r(k_x)|^2)/dk_x = 0$ можно найти проекции волнового вектора, соответствующие максимуму и минимуму отражения, окуда получаются явные выражения для коэффиента дифракционных потерь $\alpha$ и константы распространения моды $\beta$:
# \begin{equation*}
#     \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}
# \end{equation*}
# 
# \begin{equation*}
#     \beta = \frac{k_{x,max}+k_{x,min}}{2} + \frac{k_{x,max}-k_{x,min}}{2} \frac{\zeta^2-1}{\zeta^2+1}
# \end{equation*}

# #### Литература
# 
# 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)
# 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)
# 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)
# 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)
# 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)
# 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