Лицензия MIT

© Алексей Александрович Щербаков, 2024

Лекция 2.2. Фотонные кристаллы. Полный модальный базис в одномерном случае¶

В данной лекции рассмотрим аналитическое решение для одномерного фотонного кристалла, ячейка которого состоит из двух слоёв различных материалов с параметрами $\varepsilon_{1,2}$ и $\mu_{1,2}$. Данная модель является прямым аналогом модели Кронига-Пенни в электронной теории твердых тел, но в отличие от квантовомеханического случая, где периодическая система прямоугольных ям является лишь приближением потенциала, позволяющим получить аналитическое решение, в оптике кусочно-постоянное изменение материальных параметров среды описывает реально существующие структуры (если пренебрегать их внешними границами) - брэгговские зеркала.

Уравнение Гельмгольца¶

Пусть границы разделов материалов в одномерном фотонном кристалле с периодом $\Lambda$ параллельны плоскостям $XY$ декартовой системы координат, так что периодичность задается направлением оси $Z$. Без ограничения общности можно считать поля не зависящими от координаты $y$, и уравнения Максвелла разделятся для ТЕ и ТМ поляризации соответственно: \begin{equation}\tag{1} \begin{array}{c} -\dfrac{\partial E_{y}}{\partial z} = i\omega\mu H_{x} \\ \dfrac{\partial E_{y}}{\partial x} = i\omega\mu H_{z} \\ \dfrac{\partial H_{x}}{\partial z}-\dfrac{\partial H_{z}}{\partial x} = -i\omega\varepsilon E_{y} \end{array} \end{equation} и \begin{equation}\tag{2} \begin{array}{c} \dfrac{\partial H_{y}}{\partial z} = i\omega\varepsilon E_{x} \\ \dfrac{\partial H_{y}}{\partial x} = -i\omega\varepsilon E_{z} \\ \dfrac{\partial E_{x}}{\partial z}-\dfrac{\partial E_{z}}{\partial x} = i\omega\mu H_{y} \end{array} \end{equation} Уравнения Гельмгольца, следующие из этих двух, систем можно записать в едином не зависящем от поляризации виде \begin{equation}\tag{3} \dfrac{\partial}{\partial z}\left(\dfrac{1}{\eta\left(z\right)}\dfrac{\partial F}{\partial z}\right)+\dfrac{1}{\eta\left(z\right)}\dfrac{\partial^{2}F}{\partial x^{2}}+\omega^{2}\zeta\left(z\right)F=0 \end{equation} где $F^{s}\left(x,z\right) = E_{y}\left(x,z\right)$, $\eta^s=\mu$ для ТЕ поляризации и $F^p\left(x,y\right) = H_{y}\left(x,y\right)$, $\eta^p=\varepsilon$ для ТМ поляризации. Также обозначим вторую ненулевую компоненту поля в плоскости слоёв символом $G$, так что $G^s = H_x$ и $G^p = E_x$. Эту компоненту можно найти из уравнений (1),(2): $G^{x,p} = (\pm)_{s,p}(i/\omega\eta^{s,p}) dF^{s,p}/dz$. Уравнение Гельмгольца дополняется условиями непрерывности полей $F$ и $G$ на границах разделов сред между слоями кристалла и условием квазипериодичности.

Дисперсионное уравнение¶

Для решения уравнения воспользуемся методом разделения переменных. Записывая $F=u\left(x\right)v\left(z\right)$, мы приходим к двум уравнениям на новые функции: \begin{split}\tag{4} \dfrac{1}{u\left(x\right)}\dfrac{d^{2}u\left(x\right)}{dx^{2}} &= -\beta^{2} \\ \dfrac{\eta\left(z\right)}{v\left(z\right)}\dfrac{d}{dz}\left(\dfrac{1}{\eta\left(z\right)}\dfrac{\partial v\left(z\right)}{\partial z}\right) & + \left[ \omega^{2}\varepsilon\left(z\right)\mu\left(z\right)-\beta^{2} \right] = 0 \end{split} Решение первого соответствует трансляционной инвариантности полей вдоль оси $X$ (направления вдоль плоскостей кристалла) и задается экспонентой \begin{equation}\tag{5} u\left(x\right)\sim\exp\left(i\beta x\right) \end{equation} с константой распространения $\beta$. Второй оператор в случае, когда материальные параметры являются действительными функциями координат, является самосопряженным. Следовательно, его собственные значения $\beta^2$ вещественны, а собственные векторы, соответствующие этим значениям, ортогональны. Соотношение ортогональности записывается как \begin{equation}\tag{6} \dfrac{1}{\Lambda}\intop_{0}^{\Lambda}\dfrac{v_{m}\left(z\right) v^*_{n}\left(z\right)}{\omega\eta\left(z\right)}dz=\delta_{mn} \end{equation}

В рассматриваемом нами случае двух однородных слоев с толщинами $d_{1,2}$ на периоде ($d_1+d_2=\Lambda$) уравнение на функцию $v(z)$ в каждом слое сводится к уравнению второго порядка с постоянными коэффициентами \begin{equation}\tag{7} \dfrac{d^{2}v\left(z\right)}{dz^{2}}+\kappa_{1,2}^{2}v\left(z\right)=0 \end{equation} где $\kappa_{1,2}=\sqrt{\omega^{2}\varepsilon_{1,2}\mu_{1,2}-\beta^{2}}$. Решения этого уравнения выразим в форме \begin{equation}\tag{8} v\left(z\right)=\begin{cases} \mathcal{A}_{1}^{+}e^{i\kappa_{1}z}+\mathcal{A}_{1}^{-}e^{-i\kappa_{1}z} & 0\leq z\leq d_{1}\\ \mathcal{A}_{2}^{+}e^{i\kappa_{2}\left(z-d_{1}\right)}+\mathcal{A}_{2}^{-}e^{-i\kappa_{2}\left(z-d_{1}\right)} & d_{1}\leq z\leq\Lambda \end{cases} \end{equation} Постоянные коэффициенты $\mathcal{A}_{1,2}^{\pm}$ связаны друг с другом условиями сопряжения и квазипериодичности: \begin{split}\tag{9} v\left(d_{1}-0\right) = v\left(d_{1}+0\right) \\ e^{ik_{B}\Lambda}v\left(0\right) = v\left(\Lambda\right) \end{split} Те же условия, записанные для второй компоненты поля $G$, дают аналогичные условия для функции $\left(1/\eta\right)dv/dz$. В итоге получается система из четырёх уравнений \begin{split}\tag{10} \mathcal{A}_{1}^{+}e^{i\kappa_{1}d_{1}}+\mathcal{A}_{1}^{-}e^{-i\kappa_{1}d_{1}} = \mathcal{A}_{2}^{+}+\mathcal{A}_{2}^{-} \\ e^{ik_{B}\Lambda}\left(\mathcal{A}_{1}^{+} + \mathcal{A}_{1}^{-}\right) = \mathcal{A}_{2}^{+}e^{i\kappa_{2}d_{2}}+\mathcal{A}_{2}^{-}e^{-i\kappa_{2}d_{2}} \\ \dfrac{\kappa_{1}}{\omega\eta_{1}}\left(\mathcal{A}_{1}^{+}e^{i\kappa_{1}d_{1}}-\mathcal{A}_{1}^{-}e^{-i\kappa_{1}d_{1}}\right) = \dfrac{\kappa_{2}}{\omega\eta_{2}}\left(\mathcal{A}_{2}^{+}- \mathcal{A}_{2}^{-}\right) \\ e^{ik_{B}\Lambda}\dfrac{\kappa_{1}}{\omega\eta_{1}}\left(\mathcal{A}_{1}^{+}-\mathcal{A}_{1}^{-}\right) = \dfrac{\kappa_{2}}{\omega\eta_{2}}\left(\mathcal{A}_{2}^{+}e^{i\kappa_{2}d_{2}}-\mathcal{A}_{2}^{-}e^{-i\kappa_{2}d_{2}}\right) \end{split} Приравнивая детерминант к нулю, после несложных алгебраических преобразований получаем дисперсионное уравнение: \begin{equation}\tag{11} \cos\left(\kappa_{1}d_{1}\right)\cos\left(\kappa_{2}d_{2}\right)-\dfrac{1}{2}\left(\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}+\dfrac{\eta_{1}}{\kappa_{1}}\dfrac{\kappa_{2}}{\eta_{2}}\right)\sin\left(\kappa_{1}d_{1}\right)\sin\left(\kappa_{2}d_{2}\right) = \cos\left(k_{B}\Lambda\right) \end{equation} которое для действительных $\varepsilon_{1,2}$, $\mu_{1,2}$ имеет бесконечный дискретный набор решений $\beta_m^2$. Для заданных параметров кристалла имеется конечное число $\beta_m^2 > 0$, соответствующих распространяющимся вдоль оси $x$ модам и бесконечно число $\beta_m^2 < 0$, соответствующих затухающим решениям. Это можно проиллюстрировать, если построить график зависимости левой и правой частей уравнения от $\beta^2$:

Поля мод¶

Уравнения на постоянные коэффициенты $\mathcal{A}_{1,2m}^{\pm}$ можно переписать, введя новую константу: \begin{equation}\tag{12} \begin{array}{c} \mathcal{A}_{1m}^{+}=\sqrt{\mathcal{C}_{m}}\left(1-e^{ik_{B}\Lambda}e^{-i\kappa_{1m}d_{1}}e^{i\kappa_{2m}d_{2}}\right)\left(\dfrac{\eta_{2}}{\eta_{1}}\dfrac{\kappa_{1m}}{\kappa_{2m}}-1\right)\\ \mathcal{A}_{1m}^{-}=\sqrt{\mathcal{C}_{m}}\left(1-e^{ik_{B}\Lambda}e^{i\kappa_{1m}d_{1}}e^{i\kappa_{2m}d_{2}}\right)\left(\dfrac{\eta_{2}}{\eta_{1}}\dfrac{\kappa_{1m}}{\kappa_{2m}}+1\right)\\ \mathcal{A}_{2m}^{+}=\dfrac{1}{2}\sqrt{\mathcal{C}_{m}}\left[\mathcal{A}_{1m}^{+}e^{i\kappa_{1m}d_{1}}\left(1+\dfrac{\eta_{2}}{\eta_{1}}\dfrac{\kappa_{1m}}{\kappa_{2m}}\right)+\mathcal{A}_{1m}^{-}e^{-i\kappa_{1m}d_{1}}\left(1-\dfrac{\eta_{2}}{\eta_{1}}\dfrac{\kappa_{1m}}{\kappa_{2m}}\right)\right]\\ \mathcal{A}_{2m}^{-}=\dfrac{1}{2}\sqrt{\mathcal{C}_{m}}\left[\mathcal{A}_{1m}^{+}e^{i\kappa_{1m}d_{1}}\left(1-\dfrac{\eta_{2}}{\eta_{1}}\dfrac{\kappa_{1m}}{\kappa_{2m}}\right)+\mathcal{A}_{1m}^{-}e^{-i\kappa_{1m}d_{1}}\left(1+\dfrac{\eta_{2}}{\eta_{1}}\dfrac{\kappa_{1m}}{\kappa_{2m}}\right)\right] \end{array} \end{equation} Эту константу можно найти из условия ортогональности: \begin{split} 1 &= \dfrac{1}{\Lambda}\int_{0}^{\Lambda} \dfrac{v_{m}\left(z\right)v_{m}^{*}\left(z\right)}{\omega\eta\left(z\right)}dz= \\ &=\left\{ \dfrac{1}{\omega\eta_{1}}\dfrac{d_{1}}{\Lambda}\left[\left(\left|\mathcal{A}_{1m}^{+}\right|^{2}e^{-\Im\kappa_{1m}d_{1}}+\left|\mathcal{A}_{1m}^{-}\right|^{2}e^{\Im\kappa_{1m}d_{1}}\right) \dfrac{\sinh\left(\Im\kappa_{1m}d_{1}\right)}{\Im\kappa_{1m}d_{1}} + 2\Re\left(\mathcal{A}_{1m}^{+}\mathcal{A}_{1m}^{-*}e^{i\Re\kappa_{1m}d_{1}}\right)\dfrac{\sin\left(\Re\kappa_{1m}d_{1}\right)}{\Re\kappa_{1m}d_{1}}\right]+\right. \\ &\left. + \dfrac{1}{\omega\eta_{2}}\dfrac{d_{2}}{\Lambda}\left[\left(\left|\mathcal{A}_{2m}^{+}\right|^{2}e^{-\Im\kappa_{2m}d_{2}}+\left|\mathcal{A}_{2m}^{-}\right|^{2}e^{\Im\kappa_{2m}d_{2}}\right) \dfrac{\sinh\left(\Im\kappa_{2m}d_{2}\right)}{\Im\kappa_{2m}d_{2}}+2\Re\left(\mathcal{A}_{2m}^{+}\mathcal{A}_{2m}^{-*}e^{i\Re\kappa_{2m}d_{2}}\right)\dfrac{\sin\left(\Re\kappa_{2m}d_{2}\right)}{\Re\kappa_{2m}d_{2}}\right] \right\} ^{-1} \end{split}

Разложение произвольного поля по найденным собственным решениям будет иметь вид: \begin{split}\tag{13} F &= \sum_{m=1}^{\infty}v_{m}\left(z\right)\left(a_{m}^{+}e^{i\beta_{m}x}+a_{m}^{-}e^{-i\beta_{m}x}\right) \\ G=\left(\mp_{s,p}\right)\dfrac{1}{i\omega\eta}\dfrac{\partial F}{\partial x} &= \left(\mp_{s,p}\right) \sum_{m=1}^{\infty} v_{m}\left(z\right) \dfrac{\beta_{m}}{\omega\eta} \left(a_{m}^{+}e^{i\beta_{m}x}-a_{m}^{-}e^{-i\beta_{m}x}\right) \end{split} где $a_m^{\pm}$ - постоянные коэффициенты разложения.

Предел эффективной среды¶

Вернёмся к полученному дисперсионному уравнению и рассмотрим распространение волн в $\Gamma$-точке. Применяя формулы для тригонометрических функций двойного угла уравнение можно факторизовать следующим образом: \begin{split} \dfrac{\eta_{1}}{\kappa_{1}}\dfrac{\kappa_{2}}{\eta_{2}}\left[\sin\left(\dfrac{\kappa_{1}d_{1}}{2}\right)\cos\left(\dfrac{\kappa_{2}d_{2}}{2}\right)+\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}\sin\left(\dfrac{\kappa_{2}d_{2}}{2}\right)\cos\left(\dfrac{\kappa_{1}d_{1}}{2}\right)\right]\times \\ \times\left[\cos\left(\dfrac{\kappa_{1}d_{1}}{2}\right)\sin\left(\dfrac{\kappa_{2}d_{2}}{2}\right)+\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}\sin\left(\dfrac{\kappa_{1}d_{1}}{2}\right)\cos\left(\dfrac{\kappa_{2}d_{2}}{2}\right)\right]=0 \end{split} Отсюда \begin{equation*} \left[\begin{array}{c} \tan\left(\dfrac{\kappa_{1}d_{1}}{2}\right)=-\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}\tan\left(\dfrac{\kappa_{2}d_{2}}{2}\right)\\ \tan\left(\dfrac{\kappa_{2}d_{2}}{2}\right)=-\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}\tan\left(\dfrac{\kappa_{1}d_{1}}{2}\right) \end{array}\right. \end{equation*} Можно показать, что эти два типа решений соответствуют симметричным и антисимметричным модам. Переходя к пределу $k\Lambda\rightarrow 0$ подставим $\tan\alpha\approx\alpha$, что даёт: \begin{equation*} \left[\begin{array}{c} \kappa_{1}d_{1}=-\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}\kappa_{2}d_{2}\\ \kappa_{2}d_{2}=-\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}\tan\kappa_{1}d_{1} \end{array}\right.\Rightarrow\left[\begin{array}{c} \dfrac{d_{1}}{d_{2}}=-\dfrac{\eta_{2}}{\eta_{1}}\\ \kappa_{2}^{2}d_{2}\eta_{1}=-\kappa_{1}^{2}d_{1}\eta_{2} \end{array}\right. \end{equation*} Первое уравнение решений не имеет. Подставляя во второе уравнение явный вид $\kappa_{1,2}$, приходим к соотношению \begin{equation}\tag{14} \left(\omega^{2}\varepsilon_{2}\mu_{2}-\beta^{2}\right)\eta_{1}d_{2}=-\left(\omega^{2}\varepsilon_{1}\mu_{1}-\beta^{2}\right)\eta_{2}d_{1} \end{equation} В случае ТЕ поляризации для немагнитной среды получаем \begin{equation}\tag{15} \beta^{2}=\omega^{2}\mu_{0}\dfrac{\varepsilon_{1}d_{1}+\varepsilon_{2}d_{2}}{d_{1}+d_{2}} \end{equation} то есть, эффективная диэлектрическая проницаемость вычисляется как средняя проницаемость на периоде. Для ТМ поляризации \begin{equation}\tag{16} \beta^{2} = \omega^{2}\mu_{0} \left( \dfrac{\left(1/\varepsilon_{1}\right)d_{1} + \left(1/\varepsilon_{2}\right) d_{2}}{d_{2}+d_{1}} \right)^{-1} \end{equation} и вычислять нужно обратную среднюю проницаемость.

В случае $\beta=0$, когда мы хотим рассмотреть распространение волн перпендикулярно слоям кристалла, раскладываем дисперсионное уравнение по малому параметру $k_B\Lambda$, что дает не зависящее от поляризации решение \begin{equation}\tag{17} k_{B}^{2}=\omega^{2}\mu_{0}\dfrac{\varepsilon_{1}d_{1}+\varepsilon_{2}d_{2}}{d_{1}+d_{2}} \end{equation} Таким образом, фотонный кристалл в приближении эффективной среды будет вести себя как одноосный однородный анизотропный кристалл, который в главных осях описывается тензором диэлектрической проницаемости \begin{equation}\tag{18} \hat{\varepsilon}=\mathrm{diag}\left\{ \dfrac{\varepsilon_{1}d_{1}+\varepsilon_{2}d_{2}}{d_{1}+d_{2}};\thinspace\dfrac{\varepsilon_{1}d_{1}+\varepsilon_{2}d_{2}}{d_{1}+d_{2}};\thinspace\left(\dfrac{\left(1/\varepsilon_{1}\right)d_{1}+\left(1/\varepsilon_{2}\right)d_{2}}{d_{2}+d_{1}}\right)^{-1}\right\} \end{equation}

Литература¶

1. Botten, I. C., Craig, M. S., McPhedran, R. C., Adams, J. L., & Andrewartha, J. R. The Dielectric Lamellar Diffraction Grating. Optica Acta: International Journal of Optics, 28(3), 413–428 (1981)
2. Tishchenko, A.V. Phenomenological representation of deep and high contrast lamellar gratings by means of the modal method. Opt Quant Electron 37, 309–330 (2005).