Известия РАН. Серия физическая, 2022, T. 86, № 2, стр. 269-275

ВВЕДЕНИЕ

В последние годы активно исследуются многослойные метапленки и фотонные кристаллы. Эти структуры представляют значительный интерес из-за новых физических явлений, присутствующих в их оптическом отклике (см., например, [1‒3]). Предполагается использовать такие структуры в качестве элементной базы для устройств интегральной фотоники (см. например, [4–7]). В них временной масштаб оптического отклика составляет десятки фемтосекунд. Любой натурный эксперимент предваряется численным моделированием. Для этого требуются высоконадежные и эффективные численные методы.

Многие подходы (например, матричные методы, методы конечных элементов, методы дискретных источников) дают хорошие результаты для достаточно узких классов задач. Наиболее универсальны разностные методы. Для них закрепилось общее название Finite-Difference Time Domain (FDTD) method [8]. Практически все они сводятся к схеме типа “крест” [9], в которой электрическое и магнитное поля относятся к целым и полуцелым пространственным узлам и временны́м слоям. Методам FDTD посвящено более 2000 работ, но принципиальные трудности этих методов не устранены.

Во-первых, в реальных задачах всегда присутствуют границы раздела, на которых поля испытывают слабый либо сильный разрыв. Схема “крест” в этом случае включает разностное дифференцирование через границу раздела. Поэтому точность катастрофически падает, а сходимость может вовсе отсутствовать. На практике разрыв показателя преломления сглаживают. Однако это сильно ослабляет расчетное отражение.

Во-вторых, показатель преломления актуальных материалов (Si, Ge и др.) зависит от частоты сложным образом. Схемы типа FDTD позволяют учитывать только простейшие модельные законы дисперсии (полиномиальный, модель Лоренца и некоторые другие) [10]. Реальные законы дисперсии приближенно описывают комбинацией модельных законов. Это также вносит существенную погрешность [11].

В данной работе предложена новая полностью консервативная разностная схема для системы одномерных уравнений Максвелла. Она основана на принципе бикомпактности [12, 13]. Шаблон схемы содержит только два соседних пространственных узла, а границы раздела сред располагаются в узлах сетки. При этом схема явно включает условия сопряжения во всех пространственных узлах. Поэтому она реализует законы сохранения для сеточных функций сетки. Предложенная схема обеспечивает второй порядок точности даже для разрывных решений. Такие схемы для уравнений Максвелла ранее не были известны.

Для нестационарных задач предложен принципиально новый метод спектрального разложения, позволяющий учитывать дисперсию материала. Закон дисперсии может быть произвольным (в том числе заданным таблично). Этот подход имеет простую физическую интерпретацию и строгое математическое обоснование. Он применим для других линейных гиперболических задач в дисперсионных средах.

Для верификации предложенных методов проведены расчеты задачи о нормальном падении плоской волны на одномерный фотонный кристалл. Предложена специальная процедура, названная методом виртуального эксперимента, которая позволяет учитывать неидеальность изготовления образца (т.е. флуктуации толщин слоев). Этот фактор существенно влияет на интерференцию отраженных волн внутри фотонного кристалла. Рассчитан спектр прохождения через фотонный кристалл, и проведено его сравнение с ранее измеренным. Эти спектры хорошо согласуются в пределах экспериментальных погрешностей 2–5%.

КОНСЕРВАТИВНОСТЬ

Уравнения математической физики выводятся из интегральных соотношений (уравнений баланса), выражающих законы сохранения для малого объема. Разностная схема должна отражать основные свойства сплошной среды. Важнейшим свойством является выполнение разностных аналогов основных законов сохранения. Такие схемы называются консервативными [14–17]. При этом уравнение баланса для всей расчетной области должно быть алгебраическим следствием балансов в отдельных ячейках. Для построения консервативных схем в [14–16] предложен интегро-интерполяционный метод.

В сложных задачах может быть несколько законов сохранения. Примером такой является газодинамика: одновременно должны выполняться законы сохранения массы, импульса и энергии. Схемы, выражающие все нужные законы сохранения, называются полностью консервативными. Для задач электродинамики исходными физическими законами сохранения являются интегральные уравнения Максвелла. Напомним, что эти законы выражают сохранение а) потока электрической индукции через замкнутую поверхность (закон Гаусса); б) потока магнитной индукции через замкнутую поверхность (закон Гаусса для магнитного поля); в) потока магнитной индукции через незамкнутую поверхность, натянутую на контур с током (закон Фарадея); г) потока электрической индукции через незамкнутую поверхность, натянутую на замкнутый контур, по которому вычисляется циркуляция напряженности магнитного поля (теорема о циркуляции).

Дифференциальные формы основных уравнений применимы только если решение является достаточно гладким. Если же решение имеет разрыв или излом, то в точке разрыва дифференциальное уравнение неприменимо. В этом случае записывают интегральные уравнения баланса для малого объема, захватывающего области по обе стороны от разрыва. Далее этот объем стягивают в точку, совпадающую с разрывом. Это дает условия на разрыве. Например, в газодинамике таковыми являются условия Гюгонио, определяющие скачок температуры и давления в ударной волне.

В электродинамике применение интегральных уравнений Максвелла к указанному объему дает известные условия сопряжения полей на границе раздела сред. Эти условия являются прямым следствием исходных уравнений баланса. Тем самым, чтобы схема для расчета электродинамики слоистых сред была полностью консервативной, она обязательно должна учитывать не только уравнения Максвелла внутри слоев, но и условия сопряжения на границах раздела. В противном случае возникает сеточный дисбаланс, который влияет на решение во всей области.

Например, в схеме Йе метода конечных разностей во временной области (FDTD) это приводит к появлению нефизичных осцилляций и резкому падению точности (порядок точности ухудшается до первого) [18]. Для этой схемы доказано, что для однородной среды из разностных роторных уравнений следует разностный аналог дивергентных [8]. Однако эта схема не учитывает физически корректных условий сопряжения на границах раздела. Поэтому она неконсервативна согласно указанному выше определению по Тихонову-Самарскому. То же относится и к другим конечно-разностным и конечно-элементным методам, не учитывающим условия сопряжения.

СТАЦИОНАРНАЯ ЗАДАЧА

Рассмотрим распространение плоской электромагнитной волны частоты $\omega $ перпендикулярно набору из $Q$ плоско-параллельных диэлектрических пластин общей толщиной a. Обозначим координаты границ слоев через ${{\xi }_{1}} \leqslant {{\xi }_{2}} \leqslant \ldots \leqslant {{\xi }_{{Q - 1}}}.$ Пусть объемные и поверхностные токи зависят от времени как $\sim {\kern 1pt} {{e}^{{ - i\omega t}}}.$

Направим ось $z$ перпендикулярно пластинам. Тогда электрическое поле $\vec {E}$ направлено по оси x, магнитное поле $\vec {H}$ – по оси $y$, волновой вектор $\vec {k}$ – по оси z. Здесь $k = {\omega \mathord{\left/ {\vphantom {\omega c}} \right. \kern-0em} c},$ где $c$ – скорость света. Диэлектрическая проницаемость ${{\varepsilon }_{q}},$ магнитная восприимчивость ${{\mu }_{q}},$ проводимость ${{\sigma }_{q}},$ объемная плотность токов ${{\vec {J}}_{q}}$ q-й пластины, поверхностная плотность токов ${{\vec {j}}_{q}}$ на q-й границе не зависят от x и y. Вектора ${{\vec {J}}_{q}}$ и ${{\vec {j}}_{q}}$ направлены по оси x. Амплитуды полей зависят только от z. Такая задача описывается системой интегральных уравнений Максвелла

Здесь $S$ – произвольная поверхность, ограниченная контуром Γ. Используется система единиц СГС. На внешних границах поставим условия излучения Здесь ${{E}^{0}},$ ${{E}^{a}}$ – заданные амплитуды волн, падающих справа и слева. На границах слоев поставим условия сопряжения

Введем специальную сетку $\{ {{z}_{n}}\} ,$ $0 \leqslant n \leqslant N,$ $\Delta {{z}_{{n + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}} = {{z}_{{n + 1}}} - {{z}_{n}},$ у которой границы слоев являются узлами. Составим линейную интерполяцию сеточных значений полей между узлами

Здесь ${{\alpha }_{n}},$ ${{\beta }_{n}},$ ${{\gamma }_{n}},$ ${{\delta }_{n}}$ – коэффициенты интерполяции. Подставим интерполяцию (4) в задачу (1). При этом интегралы в (1) вычисляются по формуле трапеций, условия сопряжения аппроксимируются точно, а производные ${{\partial {{E}_{1}}} \mathord{\left/ {\vphantom {{\partial {{E}_{1}}} {\partial z}}} \right. \kern-0em} {\partial z}},$ ${{\partial {{E}_{Q}}} \mathord{\left/ {\vphantom {{\partial {{E}_{Q}}} {\partial z}}} \right. \kern-0em} {\partial z}}$ выразим из дифференциального уравнения Максвелла ${\text{rot}}{{E}_{q}} = i\omega {{c}^{{ - 1}}}{{\mu }_{q}}{{H}_{q}},$ записанного в узлах ${{z}_{0}}$ и ${{z}_{N}}$ соответственно. Это дает систему разностную схему для стационарной задачи

Здесь ${{\bar {z}}_{n}} = 0.5({{z}_{n}} + {{z}_{{n + 1}}}).$ Схема (5) является принципиально новой. Шаблон схемы включает только один шаг по пространству, т.е. схема является бикомпактной. Поэтому при выборе специальных сеток не возникает дифференцирования через границу раздела. Справедлива

Теорема 1. Разностная схема (5) на специальных сетках сходится со вторым порядком точности.

Доказательство. Согласно классическим теоремам теории разностных схем [16], из аппроксимации и устойчивости схемы следует ее сходимость.

Обоснуем аппроксимацию схемы (5) . Пусть внутри q-й пластины величины $\varepsilon ,$ $\mu ,$ $J$ имеют вторые непрерывные производные по координате. Тогда решение может иметь разрывы только на границах раздела, и внутри пластин имеет третью непрерывную производную. Поэтому на специальных сетках линейная интереполяция (4) аппроксимирует решение с точностью $O(\Delta {{z}^{2}}).$ Тем самым, разностная схема (5) имеет 2-й порядок аппроксимации.

Обоснуем устойчивость схемы (5) . Пусть волна заданной амплитуды распространяется в однородной среде (для определенности, в сторону увеличения $z$), и в узле ${{z}_{n}}$ ее электрическое поле равно ${{E}_{n}}.$ Можно показать, что в узле ${{z}_{{n + 1}}} = {{z}_{{n + 1}}} + \Delta {{z}_{{n + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}$ электрическое поле равно

Для реальных физических сред, у которых $\operatorname{Re} \sqrt {\varepsilon \mu } \geqslant 1,$ $\operatorname{Im} \sqrt {\varepsilon \mu } \leqslant 0$ выполняется

Если волна распространяется в слоистой среде, и в соседнем интервале $\Delta {{z}_{{n + {3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}}}}$ материальные параметры отличаются от таковых в интервале $\Delta {{z}_{{n + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}},$ то неравенство (7) тем более выполняется из-за частичного отражения волны от границы раздела сред. Аналогично можно показать, что для волны, распространяющейся в сторону уменьшения z, имеют место неравенства

Граничные условия, объемные и поверхностные токи являются источниками волн, распространяющихся в обоих направлениях оси z. Из (7) и (8) следует, что амплитуда этих волн не нарастает по мере удаления от источника. Соответственно, не будут нарастать и малые ошибки $\delta ,$ внесенные в граничные условия, объемные и поверхностные токи. В силу линейности задачи справедлива оценка

где $C$ есть мажорирующая постоянная, не зависящая от шага сетки. Отсюда следует устойчивость схемы по граничным условиям и правым частям. Теорема доказана.

Для схемы (5) применим метод сгущения сеток и оценки точности по методу Ричардсона [19–21]. В литературе предлагались различные методы локального сгущения и построения адаптивных сеток [22–24]. Однако эти подходы не дают оценок фактической точности. Напротив, глобальное сгущение сетки в методе Ричардсона дает апостериорную асимптотически точное значение погрешности. Этот подход является строго обоснованным [25].

Схема фактически построена интегро-интерполяционным методом. Поэтому для сред с непрерывно изменяющимися свойствами консервативность схемы следует из аддитивности интегралов. При наличии границ раздела консервативность обеспечивается явным учетом условий сопряжения. Таким образом, схема является полностью консервативной в смысле определения Тихонова–Самарского. Величины, сохраняемые схемой, перечислены выше.

Бикомпактность и полная консервативность обеспечивают сходимость численного решения к точному даже если последнее имеет сильные или слабые разрывы. Ранее известные схемы не обладали такими свойствами.

НЕСТАЦИОНАРНАЯ ЗАДАЧА

В постановке (1)–(3) заменим множитель $( - i\omega )$ на производную ${{\partial }_{t}}.$ Объемные ${{J}_{q}}(z,t)$ и поверхностные ${{j}_{q}}(z,t)$ токи могут произвольно зависеть от времени. Слева и справа на рассеиватель падают волновые пакеты ${{E}^{0}}(t)exp( - i\omega t),$ ${{E}^{a}}(t + {a \mathord{\left/ {\vphantom {a c}} \right. \kern-0em} c})exp( - i\omega t - ika)$ с несущей частотой $\omega $ и заданными огибающими ${{E}^{{0,a}}}.$

При распространении волнового пакета в линейной диспергирующей среде для каждой гармоники реализуются свои $\varepsilon ,$ $\mu ,$ $\sigma .$ Выполним численно преобразование Фурье падающих волновых пакетов, объемных и поверхностных токов по формуле трапеций, используя всюду одинаковые сетки по времени и наборы частот $\{ {{\omega }_{m}}\} ,$ ${{\omega }_{{m + 1}}} - {{\omega }_{m}} = \Delta {{\omega }_{m}}.$ Для каждой частоты получим свои амплитуды компонент падающих волновых пакетов $E_{m}^{{0,a}}$ и источников поля ${{({{J}_{q}})}_{m}},$ ${{({{j}_{q}})}_{m}}.$

Решим описанным выше методом набор стационарных задач, каждая из которых соответствует своей частоте ${{\omega }_{m}}.$ Полученные амплитуды просуммируем по всем частотам ${{\omega }_{m}}.$ Это и даст решение нестационарной задачи $E(z,t),$ $H(z,t).$ Это принципиально новый подход. Он был назван методом спектрального разложения. Имеет место

Теорема 2. Схема метода спектрального разложения сходится со вторым порядком точности.

Доказательство. Если функции ${{E}^{{0,a}}},$ ${{J}_{q}},$ ${{j}_{q}}$ имеют вторые непрерывные производные, то формула трапеций обеспечивает 2-й порядок аппроксимации для прямого и обратного преобразований Фурье.

Обоснуем устойчивость метода. Заметим, что в реальных задачах падающие импульсы и токи являются локализованными функциями времени. Длину промежутка времени, на котором они отличны от нуля, обозначим через T. Фурье-образы этих функций также являются локализованными. Интервал частот, на котором они отличны от нуля, обозначим через $\Omega .$

Пусть в одну из функций ${{E}^{{0,a}}},$ ${{J}_{q}},$ ${{j}_{q}}$ внесена малая ошибка $\delta .$ Тогда ошибка соответствующего фурье-образа не превосходит ${{\delta }_{\omega }} \leqslant {{(2\pi )}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}}T\max \left| \delta \right|.$ Последняя величина есть ошибка входных данных каждой из стационарных задач. Из (9) следует, что соответствующие погрешности решений стационарных задач не превосходят $\left| {\delta {{E}_{\omega }}} \right|,$ $\left| {\delta {{H}_{\omega }}} \right| \leqslant C{{(2\pi )}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}}T\max \left| \delta \right|.$ Тогда итоговая погрешность решения нестационарной задачи не превосходит $\left| {\delta E} \right|,\left| {\delta H} \right| \leqslant C{{(2\pi )}^{{ - 1}}}T\Omega \max \left| \delta \right|.$ Поэтому метод устойчив относительно малых возмущений входных данных. По классическим теоремам, отсюда следует сходимость схемы метода спектрального разложения. Теорема доказана.

Для предлагаемого метода применимы оценки погрешности по методу Ричардсона и экстраполяционное уточнение. Заметим, что сетки по z, t, $\omega $ должны сгущаться одновременно и в одинаковое число раз.

Несмотря на свою простоту, метод спектрального разложения является принципиально новым. Он позволяет учитывать произвольный закон дисперсии (в том числе заданный в виде таблицы). Этот подход рекомендуется в качестве общего метода построения разностных схем для гиперболических задач в линейной дисперсионной среде.

АПРОБАЦИЯ МЕТОДА

ЗАКЛЮЧЕНИЕ

Предложенные методы обеспечивают высокую количественную точность расчета реальных нестационарных задач фотоники. Это особенно важно при численных исследованиях сложных эффектов, таких как долгоживущие связанные состояния (поверхностные волны различных типов), а также при разработке новых устройств интегральной фотоники.

Работа выполнена при поддержке Российского научного фонда (проект № 20-71-00097).