{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9a53e7f1",
   "metadata": {},
   "source": [
    "Лицензия MIT\n",
    "\n",
    "© Алексей Александрович Щербаков, 2024"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d702e462",
   "metadata": {},
   "source": [
    "# Лекция 1.4. Резонансы в открытых системах."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d122178",
   "metadata": {},
   "source": [
    "## Резонансные состояния и квазинормальные моды\n",
    "\n",
    "В данной лекции будем рассматривать открытые оптические системы - резонаторы и волноводы.  Для краткости перепишем уравнения Максвелла в частотной области в операторной форме\n",
    "\\begin{equation}\\tag{1}\n",
    "    \\hat{\\mathcal{L}}(\\boldsymbol{r},\\omega) \\boldsymbol{F}(\\boldsymbol{r},\\omega) = \\boldsymbol{J}(\\boldsymbol{r},\\omega)\n",
    "\\end{equation}\n",
    "где $\\hat{\\mathcal{L}} = \\omega\\hat{\\mathcal{P}} - \\hat{\\mathcal{D}}$, так что\n",
    "\\begin{equation*}\n",
    "    \\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega) = \\left( \\begin{array}{cc} \\varepsilon & -i\\xi \\\\ -i\\xi^T & \\mu \\end{array} \\right), \\;\\; \\hat{\\mathcal{D}} = \\left( \\begin{array}{cc} 0 & \\nabla\\times \\\\ \\nabla\\times & 0 \\end{array} \\right)\n",
    "\\end{equation*}\n",
    "Вектор поля включает электрическую и магнитную компоненты $\\boldsymbol{F}(\\boldsymbol{r},\\omega) = (\\boldsymbol{E},\\;i\\boldsymbol{H})^T$ и составной вектор тока $\\boldsymbol{J}(\\boldsymbol{r},\\omega) = (-i\\boldsymbol{J}_e,\\;\\boldsymbol{J}_m)^T$. Тогда можно обобщить запись уравнения на функцию Грина\n",
    "\\begin{equation}\\tag{2}\n",
    "    \\hat{\\mathcal{L}}(\\boldsymbol{r},\\omega) \\mathcal{G}(\\boldsymbol{r},\\boldsymbol{r}',\\omega) = \\hat{\\mathbb{I}} \\delta(\\boldsymbol{r},\\boldsymbol{r}')\n",
    "\\end{equation}\n",
    "и интегральное решение уравнений Максвелла\n",
    "\\begin{equation}\\tag{3}\n",
    "    \\boldsymbol{F}(\\boldsymbol{r},\\omega) = \\int dV' \\mathcal{G}(\\boldsymbol{r},\\boldsymbol{r}',\\omega) \\boldsymbol{J}(\\boldsymbol{r}',\\omega)\n",
    "\\end{equation}\n",
    "\n",
    "Если возбуждение открытой системы происходит в течение короткого времени в момент $t=0$, то временной отклик системы запишется как\n",
    "\\begin{equation*}\n",
    "    \\boldsymbol{F}(\\boldsymbol{r},t) \\sim \\int d\\omega \\boldsymbol{F}(\\boldsymbol{r},\\omega) e^{-i\\omega t},\\;t>0\n",
    "\\end{equation*}\n",
    "где $\\boldsymbol{F}(\\boldsymbol{r},\\omega)$ - решение уравнения в частотной области для заданного пространственного распределения источников в нулевой момент времени. Предполагая, что анлитическое продолжение функции Грина, и, соответственно, вектора полей имеет счетное количество полюсов $\\omega_n = \\Omega_n - i\\Gamma_n$, $\\Gamma_n > 0$, основная теорема о вычетах позволяет переписать решение как\n",
    "\\begin{equation*}\n",
    "    \\boldsymbol{F}(\\boldsymbol{r},t) \\sim \\sum_n \\mathrm{Res}[\\boldsymbol{F}(\\boldsymbol{r},\\omega),\\omega_n] e^{-i\\omega_n t},\\;t>0\n",
    "\\end{equation*}\n",
    "так что затухание каждой компоненты определяется множителем $\\exp(-\\Gamma_n t)$.\n",
    "\n",
    "Поскольку правая часть уравнения на функцию Грина $(2)$ не зависит от частоты, то полюса функции Грина должны быть собственными числами уравнения\n",
    "\\begin{equation}\\tag{4}\n",
    "\\hat{\\mathcal{L}}(\\boldsymbol{r},\\omega_n) \\boldsymbol{F}_n(\\boldsymbol{r}) = 0\n",
    "\\end{equation}\n",
    "то есть, уравнения на собственные значения (моды) уравнений Максвелла в заданном открытом резонаторе или волноводе. Эти числа, вообще говоря, комплексные. Вне рассматриваемого резонатора/волновода поле представляется в виде уходящих волн, которые в трехмерном случае имеют множитель $\\exp(-i\\omega_n t + ik_nr)$ где волновое число свободного пространства, окружающего резонатор $k_n=\\omega_n\\sqrt{\\varepsilon\\mu}$. Отсюда видно, что затухающие во времени состояния открытых систем описываются экспоненциально расходящимися пространственными функциями, поэтому при разложении по таким состояниям нужно руководствоваться правилами, учитывающими данную специфику. Такие состояния в литературе называются квазинормальными модами или резонансными состояниями.\n",
    "\n",
    "Уравнение $(4)$ должно быть дополнено необходимыми граничными условиями. Решая это уравнение для комплексных частот имеют ввиду, что любые источники в заданной конфигурации рассматриваемых резонаторов локализованы, так что поле вне некоторого конечного объема имеет вид уходящих волн. Математически можно рассматривать уравнение на собственные значения как предельный случай общего уравнения $(1)$ с источником $\\boldsymbol{J}(\\boldsymbol{r},\\omega) \\sim \\omega-\\omega_n$ при $\\omega\\rightarrow\\omega_n$. При этом $\\boldsymbol{F}(\\boldsymbol{r},\\omega)\\rightarrow \\boldsymbol{F}_n(\\boldsymbol{r})$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f98f1df",
   "metadata": {},
   "source": [
    "## Разложение функции Грина по резонансным состояниям\n",
    "\n",
    "Разложение функции Грина по резнансным состояниям основано на теореме Миттаг-Леффлера: если функция $f(z)$ является аналитической всюду, кроме счетного числа полюсов $z_n$ с вычетами $R_n$, и $\\lim_{z\\rightarrow\\infty}f(z)/z^p=0$, то эта функция может быть представлена в виде\n",
    "\\begin{equation*}\n",
    "    f(z) = f_p(z) + \\sum_n\\frac{R_n}{z-z_n}\n",
    "\\end{equation*}\n",
    "где $f_p$ - многочлен степени $p-1$. Применим эту теорему для записи функции Грина, учитывая, что $\\lim_{\\omega\\rightarrow\\infty}\\mathcal{G}=0$\n",
    "\\begin{equation}\\tag{5}\n",
    "    \\mathcal{G}(\\boldsymbol{r},\\boldsymbol{r}',\\omega) = \\sum_n \\frac{\\mathcal{Q}_n(\\boldsymbol{r},\\boldsymbol{r}')}{\\omega-\\omega_n}\n",
    "\\end{equation}\n",
    "Подставим это выражение в уравнение на функцию Грина $(2)$, умножим на некотрое поле $\\boldsymbol{F}_b(\\boldsymbol{r})$, отличное от нуля в конечном объеме, и проинтегрируем обе части уравнения, что даст\n",
    "\\begin{equation*}\n",
    "    \\sum_n \\frac{\\hat{\\mathcal{L}}(\\boldsymbol{r},\\omega)\\boldsymbol{A}_n(\\boldsymbol{r})}{\\omega-\\omega_n} = \\boldsymbol{F}_b(\\boldsymbol{r})\n",
    "\\end{equation*}\n",
    "где\n",
    "\\begin{equation*}\n",
    "    \\boldsymbol{A}_n(\\boldsymbol{r}) = \\intop_{V'} dV' \\mathcal{Q}_n(\\boldsymbol{r},\\boldsymbol{r}') \\boldsymbol{F}_b(\\boldsymbol{r}')\n",
    "\\end{equation*}\n",
    "В пределе $\\omega\\rightarrow\\omega_m$ получаем систему независимых уравнений $\\hat{\\mathcal{L}}(\\boldsymbol{r},\\omega_n)\\boldsymbol{A}_n(\\boldsymbol{r}) = 0$, откуда следует, что получившиеся векторные амплитуды пропорциональны резонансным состояниям, задваемым уравнением $(4)$. А, поскольку выбор поля $\\boldsymbol{F}_b(\\boldsymbol{r})$ произволен в рамках оговоренных ограничений, то необходимо потребовать, чтобы применение оператора Грина проецировало на $n$-е резонансное состояние.\n",
    "\\begin{equation*}\n",
    "    \\mathcal{Q}_n(\\boldsymbol{r},\\boldsymbol{r}') \\sim \\boldsymbol{F}_n(\\boldsymbol{r}) \\otimes \\boldsymbol{X}(\\boldsymbol{r'})\n",
    "\\end{equation*}\n",
    "Чтобы найти вектор $\\boldsymbol{X}$, воспользуемся соотношением взаимности $\\mathcal{G}(\\boldsymbol{r}, \\boldsymbol{r}',\\omega) = \\mathcal{G}^T (\\boldsymbol{r}',\\boldsymbol{r},\\omega)$. Отсюда сразу получается, что $\\boldsymbol{X}\\sim\\boldsymbol{F}_n(\\boldsymbol{r}')$. Таким образом, выбирая соответствующую нормировку резонансных состояний, можно записать\n",
    "\\begin{equation}\\tag{6}\n",
    "    \\mathcal{G}(\\boldsymbol{r},\\boldsymbol{r}',\\omega) = \\sum_n \\frac{\\boldsymbol{F}_n(\\boldsymbol{r}) \\otimes \\boldsymbol{F}_n(\\boldsymbol{r'})}{\\omega-\\omega_n}\n",
    "\\end{equation}\n",
    "\n",
    "В ряде задач представление $(6)$ оказывается недостаточным, так что необходимо учесть разрезы на комплексной плоскости и выбор подходящих римновых поверхностей, например, за счет возникновения квадратного корня в дисперсионном уравнении для плоских волн. В более общем виде разложение тензора Грина записывается как\n",
    "\\begin{equation}\\tag{6a}\n",
    "    \\mathcal{G}(\\boldsymbol{r},\\boldsymbol{r}',\\omega) = \\sum_n \\frac{\\boldsymbol{F}_n(\\boldsymbol{r}) \\otimes \\boldsymbol{F}_n(\\boldsymbol{r'})}{\\omega-\\omega_n} + \\sum_p \\mathcal{G}_p(\\boldsymbol{r},\\boldsymbol{r}',\\omega)\n",
    "\\end{equation}\n",
    "где индекс $p$ соответствует вкладу отдельных разрезов:\n",
    "\\begin{equation*}\n",
    "    \\mathcal{G}_p(\\boldsymbol{r},\\boldsymbol{r}',\\omega) = \\frac{1}{2\\pi i} \\int_{C_p} d\\omega' \\frac{\\Delta \\mathcal{G}(\\boldsymbol{r},\\boldsymbol{r}',\\omega')}{\\omega-\\omega'}\n",
    "\\end{equation*}\n",
    "Здесь $C_m$ - путь вдоль соответствущего разреза, и $\\Delta \\mathcal{G}$ - изменение значения функции Грина на различных римановых поверхностях. Необходимо отметить, что в численных решениях вклад разрезов может быть также записан в форме $(6)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "270eb572",
   "metadata": {},
   "source": [
    "## Нормировка и полнота набора резонансных состояний\n",
    "\n",
    "Подставим полученное разложение функции Грина в объемное интегральное уравнение для вектора поля:\n",
    "\\begin{equation}\\tag{7}\n",
    "    \\boldsymbol{F}(\\boldsymbol{r},\\omega) = \\sum_n \\frac{\\boldsymbol{F}_n(\\boldsymbol{r})}{\\omega-\\omega_n}  \\int dV' \\boldsymbol{F}_n(\\boldsymbol{r'}) \\cdot \\boldsymbol{J}(\\boldsymbol{r}',\\omega)\n",
    "\\end{equation}\n",
    "Как обсуждалось выше $\\boldsymbol{F}(\\boldsymbol{r},\\omega) \\underset{\\omega\\rightarrow\\omega_{n}}{\\rightarrow}\\boldsymbol{F}_{n}(\\boldsymbol{r})$ для источников вида $\\boldsymbol{J}(\\boldsymbol{r},\\omega) = (\\omega-\\omega_n) \\tilde{\\boldsymbol{J}}(\\boldsymbol{r})$. Тогда переходя к пределу $\\omega\\rightarrow\\omega_n$ в уравнении $(7)$, получаем\n",
    "\\begin{equation}\\tag{8}\n",
    "    \\int dV' \\boldsymbol{F}_n(\\boldsymbol{r'}) \\cdot \\tilde{\\boldsymbol{J}}(\\boldsymbol{r}') = 1\n",
    "\\end{equation}\n",
    "Теперь возьмем уравнения $(1)$ и $(4)$ для поля некоторых локализованных источников и резонансных состояний, домножим первое на состояние $\\boldsymbol{F}^{\\ddagger}_{n}(\\boldsymbol{r})$, а второе на решение для данных источников $\\boldsymbol{F}(\\omega,\\boldsymbol{r})$, учтем явное представление оператора, и сложим:\n",
    "\\begin{align*}\n",
    "    \\omega \\boldsymbol{F}_{n}(\\boldsymbol{r}) \\cdot \\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega) \\boldsymbol{F}(\\boldsymbol{r},\\omega)\n",
    "    - \\omega_n \\boldsymbol{F}(\\boldsymbol{r},\\omega) \\cdot \\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega_n) \\boldsymbol{F}_{n}(\\boldsymbol{r}) + \\\\\n",
    "    + \\boldsymbol{F}(\\boldsymbol{r},\\omega) \\cdot \\hat{\\mathcal{D}}(\\boldsymbol{r},\\omega_n) \\boldsymbol{F}_{n}(\\boldsymbol{r})\n",
    "    - \\boldsymbol{F}_{n}(\\boldsymbol{r}) \\cdot \\hat{\\mathcal{D}}(\\boldsymbol{r},\\omega_n) \\boldsymbol{F}(\\boldsymbol{r},\\omega) = \\\\\n",
    "    = (\\omega-\\omega_n) \\boldsymbol{F}_{n}(\\boldsymbol{r}) \\cdot \\tilde{\\boldsymbol{J}}(\\boldsymbol{r})\n",
    "\\end{align*}\n",
    "Подставляя во вторую строчку явный вид векторов поля и оператора $\\hat{\\mathcal{D}}$, и используя разложение дивергенции векторного произведения, получаем,\n",
    "\\begin{equation*}\n",
    "    \\boldsymbol{F}(\\boldsymbol{r},\\omega) \\cdot \\hat{\\mathcal{D}}(\\boldsymbol{r},\\omega_n) \\boldsymbol{F}_{n}(\\boldsymbol{r})\n",
    "    - \\boldsymbol{F}_{n}(\\boldsymbol{r}) \\cdot \\hat{\\mathcal{D}}(\\boldsymbol{r},\\omega_n) \\boldsymbol{F}(\\boldsymbol{r},\\omega) = \n",
    "    i\\nabla \\cdot (\\boldsymbol{E}_{n} \\times \\boldsymbol{H} - \\boldsymbol{E} \\times \\boldsymbol{H}_n)\n",
    "\\end{equation*}\n",
    "Проинтегрируем общее равенство по некотрому конечному объёму, включающему локализованные источники и перейдём к пределу $\\omega\\rightarrow\\omega_n$. При этом воспользуемся разложением аналитических векторов поля и произведения $\\omega \\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega)$ продолжений вблизи точки $\\omega=\\omega_n$:\n",
    "\\begin{equation*}\n",
    "    \\omega \\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega) \\approx \\omega_n \\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega_n) + (\\omega-\\omega_n) {[\\omega \\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega)]'}_{\\omega=\\omega_n} \n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "    \\boldsymbol{E}(\\boldsymbol{r},\\omega) \\approx \\boldsymbol{E}_n(\\boldsymbol{r}) + (\\omega-\\omega_n) {[\\boldsymbol{E}(\\boldsymbol{r},\\omega)]'}_{\\omega=\\omega_n} = \\boldsymbol{E}_n(\\boldsymbol{r}) + (\\omega-\\omega_n) \\frac{1}{\\omega_n} (\\boldsymbol{r}\\nabla) \\boldsymbol{E}_n(\\boldsymbol{r})\n",
    "\\end{equation*}\n",
    "Здесь штрих обозначает производную по частоте, а последнее равенство получается, если учесть, что на границе указанного объема поле представимо в виде разложения по уходящим на бесконечность волнам и имеет место соотношение $\\boldsymbol{E}(\\boldsymbol{r},\\omega) = \\boldsymbol{E}(\\omega\\boldsymbol{r})$. В итоге учитывая соотношение $(8)$ получаем\n",
    "\\begin{equation}\\tag{9}\n",
    "    V_n + S_n = 1\n",
    "\\end{equation}\n",
    "где\n",
    "\\begin{equation}\\tag{10}\n",
    "    V_n = \\intop_V dV' \\boldsymbol{F}_n(\\boldsymbol{r}') \\cdot {[\\omega \\hat{\\mathcal{P}}(\\boldsymbol{r}',\\omega)]'}_{\\omega=\\omega_n} \\boldsymbol{F}_n(\\boldsymbol{r}') \n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\\tag{11}\n",
    "    S_n = i \\oint_{\\partial V} d\\boldsymbol{\\Sigma} (\\boldsymbol{E}_{n} \\times \\boldsymbol{H}'_n - \\boldsymbol{E}'_n \\times \\boldsymbol{H}_n)\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "На практике производная по частоте поля в резонансном состоянии может быть вычислена следующим образом. Вернемся к уравнению $(7)$ и рассмотрим его в окрестности состояния $m$, подставив соответствующее представление источника:\n",
    "\\begin{equation*}\n",
    "    \\boldsymbol{F}(\\boldsymbol{r},\\omega) = \\sum_n \\frac{\\omega-\\omega_m}{\\omega-\\omega_n} \\boldsymbol{F}_n(\\boldsymbol{r})  \\int dV' \\boldsymbol{F}_n(\\boldsymbol{r'}) \\cdot \\tilde{\\boldsymbol{J}}_m(\\boldsymbol{r}',\\omega) = \\sum_n \\frac{\\omega-\\omega_m}{\\omega-\\omega_n} \\boldsymbol{F}_n(\\boldsymbol{r})  I_{nm}\n",
    "\\end{equation*}\n",
    "Здесь $I_{mm} = 1$ согласно полученному выше соотношению. Поле вне резонансной системы раскладывается по уходящим векторным волнам. Запишем такое разложение для каждого резонансного сосотояния $\\boldsymbol{F}_n(\\boldsymbol{r}) = \\sum_L \\alpha_{Ln} \\boldsymbol{O}(\\boldsymbol{r},\\omega)$ и продифференцируем по частоте в точке $\\omega=\\omega_m$:\n",
    "\\begin{equation}\\tag{12}\n",
    "    \\boldsymbol{F}_m'(\\boldsymbol{r}) = \\sum_L [\\alpha_{Lm} \\boldsymbol{O}'(\\boldsymbol{r},\\omega_m) + \\alpha_{Lm}' \\boldsymbol{O}(\\boldsymbol{r},\\omega_m)]\n",
    "\\end{equation}\n",
    "При подстановке в соотношение ортогональности второе слагаемое не вносит вклада по следующей причине. Рассмотрим поля $\\boldsymbol{F}_{1,2}(\\boldsymbol{r},\\omega)$, возбуждаемые локализованными в области резонатора источниками $\\boldsymbol{J}_{1,2}(\\boldsymbol{r},\\omega)$. Тогда для Лоренц-взаимной среды\n",
    "\\begin{equation*}\n",
    "    \\boldsymbol{F}_1 \\cdot \\hat{\\mathcal{D}} \\boldsymbol{F}_2 - \\boldsymbol{F}_2 \\cdot \\hat{\\mathcal{D}} \\boldsymbol{F}_1 = \\boldsymbol{F}_2 \\cdot \\boldsymbol{J}_1 - \\boldsymbol{F}_1 \\cdot \\boldsymbol{J}_2\n",
    "\\end{equation*}\n",
    "Интегрирование по объему и применение теоремы Гаусса даёт\n",
    "\\begin{equation*}\n",
    "    i \\oint_{\\partial V} d\\boldsymbol{\\Sigma} (\\boldsymbol{E}_2\\times\\boldsymbol{H}_1-\\boldsymbol{E}_1\\times\\boldsymbol{H}_2) = \\int_V dV (\\boldsymbol{F}_2 \\cdot \\boldsymbol{J}_1 - \\boldsymbol{F}_1 \\cdot \\boldsymbol{J}_2)\n",
    "\\end{equation*}\n",
    "Подстановка в правую часть решений для поля через функцию Грина и применение теоремы взаимности даёт\n",
    "\\begin{equation}\\tag{13}\n",
    "    i \\oint_{\\partial V} d\\boldsymbol{\\Sigma} (\\boldsymbol{E}_2\\times\\boldsymbol{H}_1-\\boldsymbol{E}_1\\times\\boldsymbol{H}_2) = 0\n",
    "\\end{equation}\n",
    "Это уравнение показывает, что при подстановке $(12)$ в $(9)$ второе слагаемое обнуляется, и достаточно вычислять производные уходящих векторных волн по частоте, что в случае плоских, сферических или цилиндрических волн можно сделать аналитически.\n",
    "\n",
    "\n",
    "Производя вывод, такой же, как был использован для равенства $(13)$ для двух резонансных состояний $\\boldsymbol{F}_{m,n}(\\boldsymbol{r})$, получаем второе уравнение, выражающее ортогональность состояний:\n",
    "\\begin{equation}\\tag{14}\n",
    "    \\int_V dV \\boldsymbol{F}_m \\cdot [\\omega_n \\hat{\\mathcal{P}}(\\omega_n) - \\omega_m \\hat{\\mathcal{P}}(\\omega_m)] \\boldsymbol{F}_n \n",
    "    + i \\oint_{\\partial V} d\\boldsymbol{\\Sigma} (\\boldsymbol{E}_{m} \\times \\boldsymbol{H}_n - \\boldsymbol{E}_n \\times \\boldsymbol{H}_m) = 0\n",
    "\\end{equation}\n",
    "\n",
    "Переходя к вопросу о полноте системы резонансных состояний, подставим разоложение $(6)$ в уравнение для функции Грина $(2)$ и используем явный вид оператора и уравнение на состояния $(4)$. После несложных преобразований получаем\n",
    "\\begin{equation}\\tag{15}\n",
    "    \\sum_n \\frac{\\omega \\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega) - \\omega_n \\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega_n)}{\\omega-\\omega_n} \\boldsymbol{F}_n(\\boldsymbol{r}) \\otimes \\boldsymbol{F}_n(\\boldsymbol{r}') = \\hat{\\mathbb{I}} \\delta(\\boldsymbol{r}-\\boldsymbol{r}')\n",
    "\\end{equation}\n",
    "Полученное равенство выражает замыкание системы резонансных состояний и является необходимым условием полноты. В случае отсутствия дисперсии $\\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega)=\\hat{\\mathcal{P}}(\\boldsymbol{r})$ у материалов оно упрощается\n",
    "\\begin{equation}\\tag{15a}\n",
    "    \\hat{\\mathcal{P}}(\\boldsymbol{r}) \\sum_n \\boldsymbol{F}_n(\\boldsymbol{r}) \\otimes \\boldsymbol{F}_n(\\boldsymbol{r}') = \\hat{\\mathbb{I}} \\delta(\\boldsymbol{r}-\\boldsymbol{r}')\n",
    "\\end{equation}\n",
    "В более реалистичном случае дисперсии Друде-Лоренца, когда\n",
    "\\begin{equation*}\n",
    "    \\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega) = \\hat{\\mathcal{P}}_{\\infty} (\\boldsymbol{r}) + \\sum_{i} \\frac{\\hat{\\mathcal{P}}_{i}(\\boldsymbol{r})}{\\omega-\\Omega_{i}}\n",
    "\\end{equation*}\n",
    "линейная независимость лоренцианов приводит следующему выражению замыкания и дополнительным правилам сумм:\n",
    "\\begin{equation}\\tag{15b}\n",
    "    \\hat{\\mathcal{P}}_{\\infty}(\\boldsymbol{r}) \\sum_n \\boldsymbol{F}_n(\\boldsymbol{r}) \\otimes \\boldsymbol{F}_n(\\boldsymbol{r}') = \\hat{\\mathbb{I}} \\delta(\\boldsymbol{r}-\\boldsymbol{r}')\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation*}\n",
    "    \\hat{\\mathcal{P}}_{i}(\\boldsymbol{r}) \\sum_n \\frac{\\boldsymbol{F}_n(\\boldsymbol{r}) \\otimes \\boldsymbol{F}_n(\\boldsymbol{r}')}{\\omega_n-\\Omega_i} = \\hat{\\mathbb{I}} \\delta(\\boldsymbol{r}-\\boldsymbol{r}')\n",
    "\\end{equation*}\n",
    "\n",
    "Наличие такого набора соотношений может свидетельствовать об избыточности базиса, тем не менее, на практике, когда разложение производится по небольшому числу состояний, удается получить достаточную точностью расчетов. Также открытым остаётся вопрос о полноте при наличии в системе особых точек, в которых происходит одновременное вырождение и собственных частот, и собственных векторов."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c9b7c1a",
   "metadata": {},
   "source": [
    "## Поле в ближней и дальней зоне, матрица рассеяния\n",
    "\n",
    "Рассмотрим некоторый объём, поверхность которого охватыват рассматриваемый резонатор. На этой поверхности любое поле может быть разложено по полному набору приходящих и уходящих векторным волнам, $\\boldsymbol{I}_{L}(\\boldsymbol{r},\\omega)$ и $\\boldsymbol{O}_{L}(\\boldsymbol{r},\\omega)$:\n",
    "\\begin{equation}\\tag{16}\n",
    "    \\boldsymbol{F}(\\boldsymbol{r},\\omega) = \\sum_L [\\alpha^{inc}_L(\\omega) \\boldsymbol{I}_{L}(\\boldsymbol{r},\\omega) + \\alpha^{out}_L(\\omega) \\boldsymbol{O}_{L}(\\boldsymbol{r},\\omega)]\n",
    "\\end{equation}\n",
    "Набор волн может быть выбран так, чтобы удовлетворять соотношениям\n",
    "\\begin{equation*}\n",
    "     \\oint_{\\partial V} d\\boldsymbol{\\Sigma} \\cdot (\\boldsymbol{E}^{inc}_{L}\\times\\boldsymbol{H}^{out}_{L'} - \\boldsymbol{E}^{out}_{L'}\\times\\boldsymbol{H}^{inc}_{L}) = \\delta_{LL'}\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "     \\oint_{\\partial V} d\\boldsymbol{\\Sigma} \\cdot (\\boldsymbol{E}^{\\sigma}_{L}\\times\\boldsymbol{H}^{\\sigma}_{L'} - \\boldsymbol{E}^{\\sigma}_{L'}\\times\\boldsymbol{H}^{\\sigma}_{L}) = 0,\\;\\sigma=inc,out\n",
    "\\end{equation*}\n",
    "тогда коэффциенты\n",
    "\\begin{equation*}\n",
    "     \\alpha^{inc}_L = -i \\oint_{\\partial V} d\\boldsymbol{\\Sigma} \\cdot (\\boldsymbol{E}^{out}_{L}\\times\\boldsymbol{H} - \\boldsymbol{E}\\times\\boldsymbol{H}^{out}_{L})\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "     \\alpha^{sca}_L = i \\oint_{\\partial V} d\\boldsymbol{\\Sigma} \\cdot (\\boldsymbol{E}^{inc}_{L}\\times\\boldsymbol{H} - \\boldsymbol{E}\\times\\boldsymbol{H}^{inc}_{L})\n",
    "\\end{equation*}\n",
    "\n",
    "Взаимодействие резонатора с падающим полем может быть описано в терминах матрицы рассеяния:\n",
    "\\begin{equation}\\tag{17}\n",
    "     \\alpha^{out}_L= \\sum_{L'} S_{LL'} \\alpha^{inc}_{L'}\n",
    "\\end{equation}\n",
    "\n",
    "Выразим в явном виде компоненты матрицы рассеяния, выделив в физической системе некоторый фоновый оператор $\\hat{\\mathcal{P}}_b$ и добавку $\\Delta\\hat{\\mathcal{P}} = \\hat{\\mathcal{P}} - \\hat{\\mathcal{P}}_b$, описывающую локальное измерение материальных параметров за счет присутствия резонатора. Запишем полное поле в системе как сумму фонового и рассеянного $\\boldsymbol{F}_{tot} = \\boldsymbol{F}_b + \\boldsymbol{F}_{sca}$, где фоновое поле состоит из приходящих и уходящих волн и удовлетворяет уравнению $\\hat{\\mathcal{L}}_b \\boldsymbol{F}_b = 0$. Отсюда можно записать уравнение на рассеянное поле как\n",
    "\\begin{equation}\\tag{18}\n",
    "     \\hat{\\mathcal{L}}(\\boldsymbol{r},\\omega) \\boldsymbol{F}_{sca}(\\boldsymbol{r},\\omega) = -\\omega \\Delta\\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega) \\boldsymbol{F}_b\n",
    "\\end{equation}\n",
    "Решение этого уравнения через объемное интегральное уравнение с функцией Грина, разложенной по резнансным состояниям, приводит к выражению\n",
    "\\begin{equation}\\tag{19}\n",
    "     \\boldsymbol{F}_{tot}(\\boldsymbol{r},\\omega) = \\boldsymbol{F}_b(\\boldsymbol{r},\\omega) - \\sum_n \\frac{I_n(\\omega)}{\\omega-\\omega_n}\\boldsymbol{F}_n(\\boldsymbol{r})\n",
    "\\end{equation}\n",
    "где интеграл перекрытия\n",
    "\\begin{equation}\\tag{20}\n",
    "     I_n(\\omega) = \\omega \\int dV \\boldsymbol{F}_n(\\boldsymbol{r}) \\cdot \\Delta\\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega) \\boldsymbol{F}_b(\\boldsymbol{r},\\omega)\n",
    "\\end{equation}\n",
    "Это представление используется для резонансного разложения поля в ближней зоне. Для представления поля в дальней зоне перейдём к разложениям по векторным волнам: $\\boldsymbol{F}_{tot,L} = \\boldsymbol{F}_{b,L} + \\boldsymbol{F}_{sca,L}$. Соотствествующим образом можно предствить и матрицу рассеяния:\n",
    "\\begin{equation}\\tag{21}\n",
    "     S_{LL'}(\\omega) = S^b_{LL'}(\\omega) + S^{sca}_{LL'}(\\omega)\n",
    "\\end{equation}\n",
    "Используя выражения для коэффициентов приходящих и уходящих волн, можно получить\n",
    "\\begin{equation}\\tag{22}\n",
    "     S^b_{LL'}(\\omega) = i \\oint_{\\partial V} d\\boldsymbol{\\Sigma} \\cdot (\\boldsymbol{E}^{inc}_{L}\\times\\boldsymbol{H}_{b,L'} - \\boldsymbol{E}_{b,L'}\\times\\boldsymbol{H}^{inc}_{L})\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\\tag{23}\n",
    "     S^{sca}_{LL'}(\\omega) = i \\oint_{\\partial V} d\\boldsymbol{\\Sigma} \\cdot (\\boldsymbol{E}^{inc}_{L}\\times\\boldsymbol{H}_{sca,L'} - \\boldsymbol{E}_{sca,L'}\\times\\boldsymbol{H}^{inc}_{L})\n",
    "\\end{equation}\n",
    "Подставляя выражение через разность материальных операторов, приходим к выражению\n",
    "\\begin{equation}\\tag{24}\n",
    "     S_{LL'}(\\omega) = S^b_{LL'}(\\omega) - \\sum_n \\frac{\\alpha_{L}^{(n)}(\\omega)I^b_{n,L'}(\\omega)}{\\omega-\\omega_n^p}\n",
    "\\end{equation}\n",
    "где\n",
    "\\begin{equation*}\n",
    "     \\alpha^{(n)}_L = -i \\oint_{\\partial V} d\\boldsymbol{\\Sigma} \\cdot (\\boldsymbol{E}^{inc}_{L}\\times\\boldsymbol{H}_n - \\boldsymbol{E}_n\\times\\boldsymbol{H}^{inc}_{L})\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "     I^b_{n,L}(\\omega) = \\omega \\int dV \\boldsymbol{F}_n(\\boldsymbol{r}) \\cdot \\Delta\\hat{\\mathcal{P}}(\\boldsymbol{r},\\omega) \\boldsymbol{F}_{b,L}(\\boldsymbol{r},\\omega)\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20b3250c",
   "metadata": {},
   "source": [
    "## Разложение по резонансным состояниям\n",
    "\n",
    "Одним из способов нахождения разложения по резонансным состояниям сложных систем является использование известных или легко вычислимых параметров резонансных состояний в более простых системах. По аналогии с разложением задачи на фоновую и резонансную, предположим, что оператор $\\hat{\\mathcal{L}}_0$ (и его материальная часть $\\hat{\\mathcal{P}}_0$) соответствует базовой задаче с известным разложение, а её отличие от рассматриваемой задается оператором $\\delta\\hat{\\mathcal{P}}_0$. Тогда задача на собственные состояния перепишется как\n",
    "\\begin{equation*}\n",
    "    \\hat{\\mathcal{L}}_0(\\boldsymbol{r},\\omega_{\\nu}) \\boldsymbol{F}_{\\nu}(\\boldsymbol{r}) = -\\omega_{\\nu} \\delta\\hat{\\mathcal{P}}_0(\\boldsymbol{r},\\omega_{\\nu}) \\boldsymbol{F}_{\\nu}(\\boldsymbol{r})\n",
    "\\end{equation*}\n",
    "Здесь моды сложной системы нумеруются индексом $\\nu$, чтобы различать их с модами базисной простой системы. Решением уравнения будет\n",
    "\\begin{equation*}\n",
    "    \\boldsymbol{F}_{\\nu}(\\boldsymbol{r}) = -\\sum_n \\frac{\\omega_{\\nu}\\boldsymbol{F}_{n}(\\boldsymbol{r})}{\\omega_{\\nu}-\\omega_n} \\int_V dV' \\boldsymbol{F}_{n}(\\boldsymbol{r}') \\delta\\hat{\\mathcal{P}}_0(\\boldsymbol{r}',\\omega_{\\nu}) \\boldsymbol{F}_{\\nu}(\\boldsymbol{r}')\n",
    "\\end{equation*}\n",
    "Будем искать неизвестные состояния в виде разложения $\\boldsymbol{F}_{\\nu} = \\sum_n c_n \\boldsymbol{F}_{n}$. Подставляя его в обе части последнего уравнения, приходим к системе алгебраических уравнений, выражающих задачу на собственные значения:\n",
    "\\begin{equation}\\tag{25}\n",
    "    (\\omega_{\\nu}-\\omega_n) c_n = -\\omega_{\\nu} \\sum_m V_{nm}(\\omega_{\\nu}) c_m\n",
    "\\end{equation}\n",
    "где \n",
    "\\begin{equation*}\n",
    "    V_{nm} = \\int_V dV' \\boldsymbol{F}_{n}(\\boldsymbol{r}') \\delta\\hat{\\mathcal{P}}_0(\\boldsymbol{r}',\\omega_{\\nu}) \\boldsymbol{F}_{m}(\\boldsymbol{r}')\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcb72616",
   "metadata": {},
   "source": [
    "#### Литература\n",
    "\n",
    "1. S. Both and T. Weiss, [Resonant states and their role in\n",
    "nanophotonics](https://iopscience.iop.org/article/10.1088/1361-6641/ac3290), Semicond. Sci. Technol. 37 013002 (2022)\n",
    "2. Weiss, T. and Muljarov, E. A. [How to calculate the pole expansion of the optical scattering matrix from the resonant states](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.98.085433), Phys. Rev. B 98 085433 (2018)\n",
    "3. P. T. Kristensen, K. Herrmann, F. Intravaia, and K. Busch, [Modeling electromagnetic resonators using quasinormal modes](https://opg.optica.org/aop/fulltext.cfm?uri=aop-12-3-612&id=434464), Adv. Opt. Photon. 12, 612-708 (2020)"
   ]
  }
 ],
 "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
}