#
# Иллюстрация вырезанного куба из бесконечной одномерно-, двумерно- и трехмерно-периодической среды - фотонного кристалла.
#
#
#
# Терминология фотонных кристаллов возникла из аналогии с терминологией той физики твердого тела, в которой рассматриваются периодические потенциалы. Электроны в твердых кристаллических телах можно описать с помощью зонной теории. Такие понятия, как обратная решетка, зоны Бриллюэна, запрещенная зона появляются и в оптике периодических сред. Так, например, в фотонных кристаллах оказыватся возможным существование частотных интервалов, в которых отсутствуют распространяющиеся собственные решения, иначе говоря, все возбуждаемые в рамках этих интервалов волны будут затухающими. Такие интервалы называются запрещенными зонами. Возможность образования запрещенных зон для всех возможных электромагнитных волн в кристалле показывает его существенное отличие от свободного изотропного или анизотропного пространства. Поэтому можно ожидать существенного изменения характеристик различных источников электромагнитного излучения, помещенных внутри или рядом с фотонно-кристаллическими структурами. Другим направлением применения особых свойств фотонных кристаллов является создание высокодобротных резонаторов и волноводов за счет внедрения дефектов в идеальную кристаллическую решетку. При этом, если все оптические материалы, из которых состоит фотонный кристалл, являются диэлектрическими, то добротность таких резонаторов может быть очень высокой.
# ## Запрещенная зона
#
# Рассмотрим общие закономерности поведения электромагнитных волн в фотонных кристаллах. В простейшем одномерном случае немагнитных сред уравнение для собственного поля в кристалле является скалярным уравнением Гельмгольца
# \begin{equation}\tag{1}
# \left[\frac{d^2}{dz^2} + \omega^2\varepsilon(z)\mu_0\right] E(z) = 0
# \end{equation}
# с периодической функций $\varepsilon(z+n\Lambda) = \varepsilon(z)$, где $\Lambda$ - период кристалла. Согласно теореме Блоха, поле в периодическом потенциале представляется в виде произведения блоховской экспоненты на периодическую функцию. Последнюю разложим в ряд Фурье, так что
# \begin{equation}\tag{2}
# E(z) = e^{ik_0z} u(z) = e^{ik_0z} \sum_{m=-\infty}^{\infty} u_m \exp(2\pi im z/\Lambda).
# \end{equation}
# Множество точек $mK = 2\pi im/\Lambda$ задаёт сетку в обратном пространстве, где область $-K/2\leq k\leq K/2$ называется первой зоной Бриллюэна. Умножим уравнение (1) на обратную функцию диэлектрической проницаемости и также разложим эту функцию в ряд Фурье. Для простоты сначала предположим, что основной вклад в ряд дают первые три члена, так что
# \begin{equation}\tag{3}
# \varepsilon^{-1}(z) = \sum_{m=-\infty}^{\infty} \zeta_m \exp(2\pi im z/\Lambda) \approx \zeta_0 + \zeta_1 e^{2\pi i z/\Lambda} + \zeta_{-1} e^{-2\pi i z/\Lambda}.
# \end{equation}
# Подставляя (2) и (3) в (1) и используя ортогональность экспоненциальных функций, получаем систему связанных уравнений
# \begin{equation}\tag{4}
# \zeta_1 \left[ k_0 + \frac{2\pi}{\Lambda}(m-1) \right]^2 u_{m-1} + \zeta_{-1} \left[ k_0 + \frac{2\pi}{\Lambda}(m+1) \right]^2 u_{m+1} \approx \left[ \omega^2\mu_0 - \zeta_0 \left(k_0 + \frac{2\pi}{\Lambda}m\right)^2 \right] u_{m}
# \end{equation}
# Запишем уравнения для $m=0,\pm1$ и рассмотрим значения блоховского вектора на краю первой зоны Бриллюэна $k_0\approx K/2$. Если отличие диэлектрической проницаемости от свободного изотропного пространства мало, можно так же считать, что в большинстве точек k-пространства дисперсия не сильно отличается от дисперсии свободного ропстранства $\omega^2\mu_0\approx \zeta_0 k_0^2$. При этих условиях можно показать, что основной вклад в поле вносят слагаемые с $m=0,1$. Пренебрегая всеми остальными членами в (4), получием систему двух уравнений
# \begin{equation*}
# \begin{array}{c}
# \left( \omega^2\mu_0 - \zeta_0 k_0^2 \right) u_{0} - \zeta_1 \left( k_0 - \frac{2\pi}{\Lambda} \right)^2 u_{-1} \approx 0 \\
# -\zeta_{-1} k_0^2 u_{0} + \left[ \omega^2\mu_0 - \zeta_0 \left(k_0 - \frac{2\pi}{\Lambda}\right)^2 \right] u_{-1} \approx 0
# \end{array}
# \end{equation*}
# Обозначим $\Delta k = k-\pi/\Lambda$, так что $|\Delta k|\ll \pi/\Lambda$. Приравнивая детерминант системы к нулю и решая относительно собственной частоты, находим два решения
# \begin{equation*}
# \omega_{\pm} \approx \frac{K}{2} \sqrt{\zeta_0\pm|\zeta_1|} \pm \frac{1}{\sqrt{\zeta_0}|\zeta_1|} \frac{(\Delta k)^2}{K/2} \left( \zeta_0^2 -\frac{\zeta_1^2}{4} \right)
# \end{equation*}
# Запрещенная зона существует в области частот $(K/2)\sqrt{\zeta_0-|\zeta_1|} < \omega < (K/2)\sqrt{\zeta_0+|\zeta_1|}$.
#
# При отсутствии модуляции $\zeta_m=0$ для всех $m\neq0$ и запрещенная зона исчезает - получается так называемое приближение пустой решетки. При ненулевой модуляции существует бесконечное число запрещенных зон - выше приведена оценка для самой нижней. Пары совпадающих решений, соответствующие пересечениям дисперсионных линий при пустой решетке, расщепляются. Поскольку при наличии периодичности решения, отличающиеся на $K=2\pi/\Lambda$ эквивалентны, решения на краях зон Бриллюэна соответствуют интерференции волн, распространяющихся в противоположных направлениях, что приводит к образованию стоячих волн, групповая скорость которых $v_g = d\omega/dk_0$ равна нулю, а, значит, касательная к дисперсионной кривой $\omega(k_0)$ на границах имеет нулевой угол наклона.
# In[4]:
# Первая зона зона Бриллюэна одномерного фотонного кристалла и приближение пустой решетки
bd1.plot_1D_PhC_dispersion(eps1=1, eps2=2.5, d1=0.1, d2=0.2, n=5, nb=4)
# ## Аномалия групповой скорости
#
# В случае двумерных кристаллов существует два неколлинеарных направления периодичности и, соответственно, два базисных вектора обратной решетки. Первая зона Бриллюэна является ячейкой двумерной Вигнера-Зейца в обратном пространстве. Точки высокой симметрии этой зоны обычно обозначаются заглавными греческими буквами. Центр первый зоны Брилююэна называется $\Gamma$-точкой. Дисперсионные диаграммы изображают с помощью линий, соответствующих прямым траекториям, соединяющим точки высокой симметрии.
#
#
#
#
# Первая зона зона Бриллюэна двумерных кристаллов с квадратной и гексагональной решеткой.
#
#
#
# В приближении пустой решетки дисперсия зон задаётся уравнением
# \begin{equation}\tag{5}
# \left(k_x + mK_x\right)^2 + \left(k_y + nK_y\right)^2 = \omega^2\varepsilon\mu_0
# \end{equation}
# Например, для квадратной решетке два периода задаются векторами $\boldsymbol{p}_1 = (\Lambda,0)^T$ и $\boldsymbol{p}_2 = (0,\Lambda)^T$. Соответствующие им векторы обратной решетки $\boldsymbol{K}_1 = (2\pi/\Lambda,0)^T$ и $\boldsymbol{K}_2 = (0,2\pi/\Lambda)^T$. Траектория $\Gamma X$ задаётся изменением вектора $\boldsymbol{k} = (k_x,0)^T$. Для третьей зоны с $m=2$ можно видеть, что вдоль всей траектории наклон дисперсионной кривой остается весьма малым, что соответствует малой групповой скорости волн. Этот эффект называется аномалией групповых скоростей. В научных публикациях были продемонстрированы структуры фотонных кристаллов, допускающих появления ветвей дисперсионных диаграмм с почти нулевой групповой скоростью.
# In[1]:
from pathlib import Path
import sys
import os
path_root = os.path.dirname(os.path.dirname(os.path.abspath('')))
sys.path.append(str(path_root))
import ANMOP.code.PhC2D_empty_lattice as ela
# In[2]:
ela.plot_empty_lattice_dispersion()
# ## Скейлинг и симметрия относительно обращения времени
#
# Если нормировать волновые векторы на модуль вектора обратной решетки $\boldsymbol{k}\rightarrow (1/K) \boldsymbol{k} = (\Lambda/2\pi) \boldsymbol{k}$, а частоту измерять в единицах $2\pi c/\Lambda$, то при условии незменности значений материальных параметров электромагнитные свойства фотонных кристаллов как решений уравнений Максвелла не изменятся. Это следует из подстановки в уравнения координт и времени, нормированных на период и $\Lambda/c$ соответственно. Это свойство называется скейлингом.
#
# Волновое уравнение не зависит от замены знака времени $t'=-t$, поэтому решение $\tilde{\boldsymbol{E}}(\boldsymbol{r},t') = \boldsymbol{E}(\boldsymbol{r},-t)$ также удовлетворяет этому уравнению. Следовательно, собственные частоты $\omega_{\boldsymbol{k}n}$ - решения задачи на собственные значения - остаются прежниними, поэтому
# \begin{equation*}
# \tilde{\boldsymbol{E}}(\boldsymbol{r},t') = \boldsymbol{u}_{\boldsymbol{k}n}(\boldsymbol{r}) \exp[i(\boldsymbol{k}\boldsymbol{r} - \omega_{\boldsymbol{k}n}t')]
# \end{equation*}
# Обратное преобразование даёт
# \begin{equation*}
# \boldsymbol{E}(\boldsymbol{r},t) = \left( \boldsymbol{u}^*_{\boldsymbol{k}n}(\boldsymbol{r}) \exp[i(-\boldsymbol{k}\boldsymbol{r} - \omega_{\boldsymbol{k}n}t)] \right)^*
# \end{equation*}
# Поскольку измеряемое поле является действительной частью комплексного решения, можно видеть, что
# \begin{equation}\tag{6}
# \begin{array}{c}
# \omega_{-\boldsymbol{k}n} = \omega_{\boldsymbol{k}n}, \\
# \boldsymbol{u}_{-\boldsymbol{k}n} = \boldsymbol{u}^*_{\boldsymbol{k}n}
# \end{array}
# \end{equation}
# Отсюда следует, в частности, симметрия дисперсионных диаграмм относительно инверсии в сопряженном пространстве.
# ### Литература
#
# 1. K. Sakoda, [Optical Properties of Photonic Crystals](https://materias.df.uba.ar/introNanoFot2014/files/2012/07/Optical-Properties-of-Photonic-Crystals_sakoda.pdf), Springer Science & Business Media (2004)
# 2. J.D. Joannopoulos, S.G. Johnson, J.N. Winn, and R.D. Meade, [Photonic Crystals. Molding the Flow of Light (Second Edition)](), Princeton University Press (2011)
# 
#