Символьно-численная реализация модели адиабатических волноводных мод для двумерных нерегулярных волноводов

Полный текст

1. Введение

В работе строятся символьно-численные решения уравнений Максвелла, описывающие волноводные моды двумерного плавно-нерегулярного волновода в рамках модели адиабатических волноводных мод (АВМ) [1, 2, 3]. В основе модели АВМ лежит приближение коротких волн [4], адаптированное для волноводного распространения, которое представляет собой асимптотический ряд по малой величине, представляющей собой обратную частоту. Асимптотические методы удобны тем, что нулевое приближение отыскивается, как правило, в символьном виде. Учитывая асимптотический характер решения, в модели АВМ удается получить ряд промежуточных результатов в символьном виде – поэтому компьютерная алгебра используется в качестве основного инструмента исследования.

В нулевом приближении асимптотического разложения модели АВМ уравнения Максвелла символьно редуцируются к системе четырех обыкновенных дифференциальных уравнений первого порядка и двум дополнительным соотношениям. Систему дифференциальных уравнений удается решить символьно [5] в каждом слое многослойного волновода. Используя символьные представления решений формулируется система граничных уравнений в виде однородной системы линейных алгебраических уравнений (СЛАУ) с символьной матрицей коэффициентов. Однородная СЛАУ с символьной матрицей коэффициентов решается символьно в настоящей работе. Условие разрешимости однородной СЛАУ формулируется в виде нелинейного уравнения, которое решается в работе численно методом продолжения по параметру.

2. Методы

2.1. Модель адиабатических волноводных мод в нулевом приближении

В работе [5] получен нулевой вклад в адиабатическое приближение волноводного решения уравнений Максвелла вида:

$\left\{\begin{array}{c}\overrightarrow{E}\left(x,z,t\right)\\ \overrightarrow{H}\left(x,z,t\right)\end{array}\right\}=\left\{\begin{array}{c}{\overrightarrow{E}}_0\left(x;z\right)\\ {\overrightarrow{H}}_0\left(x;z\right)\end{array}\right\}=\exp \left\{i\omega t-i{k}_0\phi \left(z\right)\right\},$ (1)

причем

$\left\{\begin{array}{c}\epsilon \frac{\partial E_0^y}{\partial x}=-i{k}_0\epsilon \mu H_0^z,\\ \epsilon \frac{\partial E_0^z}{\partial x}=i{k}_0\left(\epsilon \mu -{\left(\phi \text{'}\left(z\right)\right)}^2\right)H_0^y,\\ \mu \frac{\partial H_0^y}{\partial x}=i{k}_0\epsilon \mu E_0^z,\\ \mu \frac{\partial H_0^z}{\partial x}=-i{k}_0\left(\epsilon \mu -{\left(\phi \text{'}\left(z\right)\right)}^2\right)E_0^y,\end{array}\right.$ (2)

$E_0^x=\frac{1}{\epsilon }\phi \text{'}\left(z\right)H_0^y,H_0^x=-\frac{1}{\mu }\phi \text{'}\left(z\right)E_0^y.$  (3)

Систему уравнений (2)–(3) следует дополнить условиями сопряжения электромагнитного поля на границах раздела сред [6] для рассматриваемого многослойного волновода: на границах раздела диэлектрических сред выполняются граничные условия сопряжения полей

${\left[n\times \overrightarrow{E}\right]}_{x=\mathrm{0,}h\left(z\right)}=\overrightarrow{0},{\left[n\times \overrightarrow{H}\right]}_{x=\mathrm{0,}h\left(z\right)}=\overrightarrow{0},$ (4)

где через ${[ \overrightarrow{f} ]}_{x = 0,h\left(z\right)}$ обозначен скачок векторной величины $\overrightarrow{f}$ на границах x = 0,h(z). Кроме того, выполняются асимптотические граничные условия на бесконечности [6]:

$||\overrightarrow{E}||\underset{x\to \pm \infty }{\to }\mathrm{0,}||\overrightarrow{H}||\underset{x\to \pm \infty }{\to }0.$ (5)

2.2. Геометрия рассматриваемой структуры

Рассматривается трехслойный плавно-нерегулярный волновод, геометрия которого представлена на рис. 1.

 

Рис. 1. Геометрия двумерного плавно-нерегулярного волноводного перехода между двумя регулярными волноводами.

 

Параметры волновода следующие: μ_c = μ_f = μ_s = 1, ε_c = 1, ε_f = 1.565², ε_s = 1.47², толщины определены как h₁ = 2λ, h₂ = 3λ, а L = 100λ, где λ – длина волны, λ = 0.55 [мкм]. Переменная толщина h(z) определена следующим образом:

$h\left(z\right)=2\left(h_1-h_2\right){\left(\frac{z}{L}\right)}^3-3\left(h_1-h_2\right){\left(\frac{z}{L}\right)}^2+h_1,$ (6)

причем для h′(z) выполняется |h′(z)| << 1, то есть h′(z) является малым параметром при каждом фиксированном z.

В работе вычисляются адиабатические волноводные моды для описанной выше структуры в символьно-численном виде.

2.3. Символьный метод решения задачи

В работе система (2) решается в символьном виде. Коэффициенты системы (2) заданы в символьном виде, причем ε, μ для рассматриваемого многослойного волновода есть кусочно-постоянные функции – известные константы для каждого из слоев. В каждом слое систему (2) решаем символьно в системе компьютерной алгебры Maple с помощью функции dsolve [7] и получаем разложение решения по фундаментальной системе решений (ФСР) с неопределенными коэффициентами. Решения в полубесконечных слоях удовлетворяют условиям (5), поэтому константы, стоящие перед нарастающими на бесконечности функциями ФСР, будут определены и равны нулю.

Условия непрерывности (4) записываем символьно в системе компьютерной алгебры Maple [7] и получаем систему граничных уравнений вида

$M\left(z,\gamma \left(z\right)\right)\overrightarrow{C}\left(z\right)=\overrightarrow{0},$ (7)

где γ(z) = j′(z), вектор $\overrightarrow{C}$(z) определяет константы разложения решения по ФСР в каждом слое при фиксированном z. Условие разрешимости системы (7) есть равенство нулю определителя системы

$\det M\left(z,\gamma \left(z\right)\right)=0.$ (8)

Система линейных алгебраических уравнений (7) с символьной матрицей коэффициентов решается в работе символьно методом из раздела 3.3 работы [8].

Cимвольное выражение детерминанта определяется с помощью функции Determinant пакета LinearAlgebra. Уравнение (8) характеризуется вещественными корнями при z ≤ 0 и при z ≥ L, поэтому корни для z = 0 отыскиваются с помощью встроенной в Maple функции RootOf [7]. В интервале 0 < z < L корни уравнения (8) могут быть комплексными, поэтому для их отыскания в системе компьютерной алгебры использован метод продолжения по параметру, описанный ниже.

Описанный метод решения нелинейного уравнения (8), а также символьное решение системы (7) реализованы в системе компьютерной алгебры Maple [7].

2.4. Метод продолжения по параметру

Рассмотрим уравнение F(x,y) = 0 относительно искомой величины y, которое содержит параметр x. Пусть это уравнение имеет при a ≤ x ≤ b решение y(x) и при x = a такое решение y(a) = y_a известно, то есть F(a,y)_a = 0. Тогда в пространстве (x, y) точка y, соответствующая решениям рассматриваемого уравнения, описывает непрерывную кривую k, проходящую через точки (a, y_a) и (b, y_b), где F(b, y_b) = 0.

Идея метода продолжения по параметру [9] состоит в построении решения отталкиваясь от (a, y_a) и двигаясь вдоль кривой k.

Дифференцирование рассматриваемого уравнения по параметру приводит к дифференциальному уравнению

$\frac{dy}{dx}=-\frac{\partial F/\partial x}{\partial F/\partial y},$ (9)

которое вместе с условием

$y\left(a\right)=y_a$ (10)

образует задачу Коши.

В работе дифференциальное уравнение (9) формулируется в системе компьютерной алгебры с помощью встроенной в Maple функции diff. Начальное условие (10) формулируется численно в системе компьютерной алгебры с помощью непостредственного решения уравнения встроенной функцией RootOf.

Сформулированная таким образом задача Коши решается численно в системе компьютерной алгебры Maple с помощью функции dsolve методом rkf45 с разными abserr, relerr.

3. Результаты

3.1. Символьное исследование системы

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

$M=\left(\begin{array}{cc}{M}_1& 0\\ 0& M_2\end{array}\right),\overrightarrow{C}=\left(\begin{array}{c}{\overrightarrow{C}}_1\\ {\overrightarrow{C}}_2\end{array}\right).$ (11)

и в результате система (7) разбиваеся на две независимых системы M₁ $\overrightarrow{C}$₁ = $\overrightarrow{0}$ и M₂ $\overrightarrow{C}$₂ = $\overrightarrow{0}$.

Матрица M₁ имеет следующую структуру

$M_1=\left(\begin{array}{cccc}{m}_{\mathrm{1,1}}^1& m_{\mathrm{1,2}}^1& m_{\mathrm{1,3}}^1& 0\\ 0& m_{\mathrm{2,2}}^1& m_{\mathrm{2,3}}^1& m_{\mathrm{2,4}}^1\\ 0& m_{\mathrm{3,2}}^1& m_{\mathrm{3,3}}^1& m_{\mathrm{3,4}}^1\\ m_{\mathrm{4,1}}^1& m_{\mathrm{4,2}}^1& m_{\mathrm{4,3}}^1& 0\end{array}\right),$  (12)

элементы матрицы определены следующим образом

$m_{\mathrm{1,1}}^1=ih\text{'}{\eta }_c\phi \text{'}-{\eta }_c^2,$ (13)

$m_{\mathrm{1,2}}^1=ih\text{'}{\eta }_f\phi \text{'}{e}^{{\eta }_f{k}_{0}h}-{\eta }_f^2{e}^{{\eta }_f{k}_{0}h},$ (14)

$m_{\mathrm{1,3}}^1=-ih\text{'}{\eta }_f\phi \text{'}{e}^{-{\eta }_f{k}_{0}h}+{\eta }_f^2{e}^{-{\eta }_f{k}_{0}h},$ (15)

$m_{\mathrm{2,2}}^1=i{\eta }_f{\epsilon }_f,$ (16)

$m_{\mathrm{2,3}}^1=-i{\eta }_f{\epsilon }_f,$ (17)

$m_{\mathrm{2,4}}^1=-i{\eta }_s{\epsilon }_s,$ (18)

$m_{\mathrm{3,2}}^1=-{\eta }_f^2,$ (19)

$m_{\mathrm{3,3}}^1=-{\eta }_f^2,$ (20)

$m_{\mathrm{3,4}}^1={\eta }_s^2,$ (21)

$m_{\mathrm{4,1}}^1=-i{\eta }_c{\epsilon }_c,$ (22)

$m_{\mathrm{4,2}}^1=-i{\eta }_f{\epsilon }_f{e}^{{\eta }_f{k}_{0}h},$ (23)

$m_{\mathrm{4,3}}^1=i{\eta }_f{\epsilon }_f{e}^{-{\eta }_f{k}_{0}h}.$ (24)

Матрица M₂ определена ниже

$M_2=\left(\begin{array}{cccc}0& m_{\mathrm{1,2}}^2& m_{\mathrm{1,3}}^2& m_{\mathrm{1,4}}^2\\ 0& m_{\mathrm{2,2}}^2& m_{\mathrm{2,3}}^2& m_{\mathrm{2,4}}^2\\ m_{\mathrm{3,1}}^2& m_{\mathrm{3,2}}^2& m_{\mathrm{3,3}}^2& 0\\ m_{\mathrm{4,1}}^2& m_{\mathrm{4,2}}^2& m_{\mathrm{4,3}}^2& 0\end{array}\right),$  (25)

элементы матрицы определены следующим образом

$m_{\mathrm{1,2}}^2=-{\eta }_f^2,$ (26)

$m_{\mathrm{1,3}}^2=-{\eta }_f^2,$ (27)

$m_{\mathrm{1,4}}^2={\eta }_s^2,$ (28)

$m_{\mathrm{2,2}}^2=i{\eta }_f{\mu }_f,$ (29)

$m_{\mathrm{2,3}}^2=-i{\eta }_f{\mu }_f,$ (30)

$m_{\mathrm{2,4}}^2=i{\eta }_s{\mu }_s,$ (31)

$m_{\mathrm{3,1}}^2=ih\text{'}{\eta }_c\phi \text{'}-{\eta }_c^2,$ (32)

$m_{\mathrm{3,2}}^2=-ih\text{'}{\eta }_f\phi \text{'}{e}^{-{\eta }_f{k}_{0}h}+{\eta }_f^2{e}^{-{\eta }_f{k}_{0}h},$ (33)

$m_{\mathrm{3,3}}^2=ih\text{'}{\eta }_f\phi \text{'}{e}^{{\eta }_f{k}_{0}h}+{\eta }_f^2{e}^{{\eta }_f{k}_{0}h},$ (34)

$m_{\mathrm{4,1}}^2=i{\eta }_c{\mu }_c,$ (35)

$m_{\mathrm{4,2}}^2=-i{\eta }_f{\mu }_f{e}^{-{\eta }_f{k}_{0}h},$ (36)

$m_{\mathrm{4,3}}^2=i{\eta }_f{\mu }_f{e}^{{\eta }_f{k}_{0}h},$ (37)

где ${\eta }_{\alpha }=\sqrt{{{\phi }^{\text{'}}}^2-{\epsilon }_{\alpha }{\mu }_{\alpha }}$, a = c, f, s. Решение подсистемы M₁ $\overrightarrow{C}$₁ = $\overrightarrow{0}$ получаем символьно, используя метод из раздела 3.3 работы [8]:

$\begin{array}{c}{C}_{\mathrm{1,1}}=\frac{{\eta }_f^2{\eta }_s\left({\eta }_s{\epsilon }_f-{\eta }_f{\epsilon }_s\right)\left(i{\eta }_f+h\text{'}\phi \text{'}\right)e^{-{\eta }_f{k}_{0}h}}{{\eta }_c\left(i{h}^{\text{'}}{\phi }^{\text{'}}-{\eta }_c\right)}+\\ +\frac{{\eta }_f^2{\eta }_s\left({\eta }_s{\epsilon }_f+{\eta }_f{\epsilon }_s\right)\left(i{\eta }_f-h^{\text{'}}{\phi }^{\text{'}}\right)e^{{\eta }_f{k}_{0}h}}{{\eta }_c\left(i{h}^{\text{'}}{\phi }^{\text{'}}-{\eta }_c\right)},\end{array}$ (38)

$C_{\mathrm{1,2}}=-i{\eta }_f{\eta }_s\left({\eta }_s{\epsilon }_f+{\eta }_f{\epsilon }_s\right),$ (39)

$C_{\mathrm{1,3}}=-i{\eta }_f{\eta }_s\left({\eta }_s{\epsilon }_f-{\eta }_f{\epsilon }_s\right),$ (40)

$C_{\mathrm{1,4}}=-2i{\eta }_f^3{\epsilon }_f.$ (41)

Решение подсистемы M₂ $\overrightarrow{C}$₂ = $\overrightarrow{0}$ получаем аналогичным образом:

$\begin{array}{c}{C}_{\mathrm{2,1}}=-\frac{i{\eta }_f^2{\eta }_s{\mu }_f\left({\eta }_f{\mu }_s-{\eta }_s{\mu }_f\right)e^{-{\eta }_f{k}_{0}h}}{{\eta }_c{\mu }_c}-\\ -\frac{i{\eta }_f^2{\eta }_s{\mu }_f\left({\eta }_f{\mu }_s+{\eta }_s{\mu }_f\right)e^{{\eta }_f{k}_{0}h}}{{\eta }_c{\mu }_c},\end{array}$ (42)

$C_{\mathrm{2,2}}=i{\eta }_f{\eta }_s\left({\eta }_s{\mu }_f-{\eta }_f{\mu }_s\right),$ (43)

$C_{\mathrm{2,3}}=i{\eta }_f{\eta }_s\left({\eta }_s{\mu }_f+{\eta }_f{\mu }_s\right),$ (44)

$C_{\mathrm{2,4}}=2i{\eta }_f^3{\mu }_f.$ (45)

На основе полученных символьных решений (38)–(45) определяем символьное представление компонент электромагнитного поля в каждом из трех слоев рассматриваемого волновода.

Амплитудная часть электромагнитного поля в покровном слое (x > h(z)) определяется следующим образом:

${\overrightarrow{E}}_0^c=\left(\begin{array}{c}-i{\eta }_c\phi \text{'}{C}_{\mathrm{1,1}}\\ i{\eta }_c{\mu }_c{C}_{\mathrm{2,1}}\\ {\eta }_c^2{C}_{\mathrm{1,1}}\end{array}\right)e^{-{\eta }_c{k}_0\left(x-h\right)},$ (46)

${\overrightarrow{H}}_0^c=\left(\begin{array}{c}-i{\eta }_c\phi \text{'}{C}_{\mathrm{2,1}}\\ -i{\eta }_c{\epsilon }_c{C}_{\mathrm{1,1}}\\ {\eta }_c^2{C}_{\mathrm{2,1}}\end{array}\right)e^{-{\eta }_c{k}_0\left(x-h\right)}.$ (47)

Амплитудная часть электромагнитного поля в волноводном слое (0 < x < h(z)) определяется следующим образом:

${\overrightarrow{E}}_0^f=\left(\begin{array}{c}i{\eta }_f\phi \text{'}{C}_{\mathrm{1,2}}\\ -i{\eta }_f{\mu }_f{C}_{\mathrm{2,3}}\\ {\eta }_f^2{C}_{\mathrm{1,2}}\end{array}\right)e^{{\eta }_f{k}_{0}x}\text{}+\left(\begin{array}{c}-i{\eta }_f\phi \text{'}{C}_{\mathrm{1,3}}\\ i{\eta }_f{\mu }_f{C}_{\mathrm{2,2}}\\ {\eta }_f^2{C}_{\mathrm{1,3}}\end{array}\right)e^{-{\eta }_f{k}_{0}x},$ (48)

${\overrightarrow{H}}_0^f=\left(\begin{array}{c}i{\eta }_f\phi \text{'}{C}_{\mathrm{2,3}}\\ i{\eta }_f{\epsilon }_f{C}_{\mathrm{1,2}}\\ {\eta }_f^2{C}_{\mathrm{2,3}}\end{array}\right)e^{{\eta }_f{k}_{0}x}+\left(\begin{array}{c}-i{\eta }_f\phi \text{'}{C}_{\mathrm{2,2}}\\ -i{\eta }_f{\epsilon }_f{C}_{\mathrm{1,3}}\\ {\eta }_f^2{C}_{\mathrm{2,2}}\end{array}\right)e^{-{\eta }_f{k}_{0}x}.$ (49)

Амплитудная часть электромагнитного поля в подложке (x < 0) определяется следующим образом:

${\overrightarrow{E}}_0^s=\left(\begin{array}{c}i{\eta }_s\phi \text{'}{C}_{\mathrm{1,4}}\\ -i{\eta }_s{\mu }_s{C}_{\mathrm{2,4}}\\ {\eta }_s^2{C}_{\mathrm{1,4}}\end{array}\right)e^{{\eta }_s{k}_{0}x},$ (50)

${\overrightarrow{H}}_0^s=\left(\begin{array}{c}i{\eta }_s\phi \text{'}{C}_{\mathrm{2,4}}\\ i{\eta }_s{\epsilon }_s{C}_{\mathrm{1,4}}\\ {\eta }_s^2{C}_{\mathrm{2,4}}\end{array}\right)e^{{\eta }_s{k}_{0}x}.$ (51)

В приведенных символьных выражениях (46)–(51) только для величины φ′ неизвестно символьное представление, так как она является решением нелинейного уравнения (8) и определяется численно.

3.2. Численное решение нелинейного уравнения

Используем для решения нелинейного уравнения (8) метод продолжения по параметру, описанный в разделе 1 настоящей работы. Вычисленные величины ϒ_j(z), j = 1..4 приведены на рис. 2, 3. Приведем также невязки решений, полученных методом продолжения по параметру. Величины ${\delta }_j\left(z\right)=\left|\det M(z,{\gamma }_j(z\left)\right)\right|$, j = 1..4 (невязки) приведены на рис. 4, 5.

 

Рис. 2. Величины ℜ(γ_j(z)), j = 1..4 для z ∈ [0,L].

 

Рис. 3. Величины ℑ(γ_j(z)), j = 1..4 для z ∈ [0,L].

 

Рис. 4. Невязки δ_j(z), j = 1..4 для z ∈ [0,L] при значениях abserr и relerr, равных 10^–10.

 

Рис. 5. Невязки δ_j(z), j = 1..4 для z ∈ [0,L] при значениях abserr и relerr, равных 10^–12.

 

4. Обсуждение

В работе получены символьные представления направляемых мод плавно-нерегулярного волноводного перехода (рис. 1) между двумя регулярными волноводами в рамках модели адиабатических волноводных мод. Символьные представления направляемых мод во всех слоях волновода достаточно компактны (46)–(51) и позволяют построить волноводные моды в символьно-численном виде, где численно отыскивается только ϒ(z) = φ′(z). Кроме того, символьный вид направляемых мод модели АВМ в нулевом приближении асимптотического разложения (46)–(51) будет использован при построении первого приближения асимптотического метода.

Важно отметить, что полученные символьные выражения (46)–(51) не единственны. Каждая подсистема M₁ $\overrightarrow{C}$₁ = $\overrightarrow{0}$ и M₂ $\overrightarrow{C}$₂ = $\overrightarrow{0}$ имеет по 4 различных символьных представлений решения, как в работе [8]. Вектор решения $\overrightarrow{C}$, определенный в (11), будет иметь 16 различных символьных представлений. Кроме того, любая линейная комбинация этих 16 символьных представлений также будет решением однородной системы (7). Вопрос выбора наиболее удобного символьного представления решения остается открытым.

В рамках приближенного вычисления корней нелинейного уравнения вычислены функции γ_j(z), описывающие переменные коэффициенты фазового замедления для j-й волноводной моды. В работе [10] приближенно вычислялись функции γ_j(z) методом разложения по малому параметру h′(z), отвечающему малому наклону криволинейной границы раздела. В работе [10] построены численно нулевое, первое и второе приближения по малому параметру и вычислены невязки. В обеих работах рассматривается одинаковая структура.

В настоящей работе вместо разложения по малому параметру использован метод продолжения по параметру, который дает невязку, не превышающую 4.5 × 10^–11, согласно рис. 4, и невязку, не превышающую 4.5 × 10^–13 согласно рис. 5. В работе [10] получена меньшая точность порядка 10^–8. Другими словами, точность метода продолжения по параметру на несколько порядков выше, чем у метода малого параметра из [10], и меняется в зависимости от значений параметров abserr, relerr метода rkf45. В методе малого параметра из [10] решение представляется в виде разложения, коэффициенты которого известны в виде довольно громоздких символьных выражений. Обобщая вышесказанное, метод малого параметра менее удобен для решения нелинейного уравнения (8), чем использованный в настоящей работе метод продолжения по параметру.

5. Заключение

Модель адиабатических волноводных мод позволяет сформулировать задачу расчета направляемых мод в символьном виде и решить ее в символьно-численном виде.

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

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

Источник финансирования

Публикация выполнена при поддержке Программы стратегического академического лидерства РУДН, проект № 021934-0-000.

Об авторах

Д. В. Диваков

Российский университет дружбы народов

Автор, ответственный за переписку.
Email: divakov_dv@pfur.ru
Россия, 117198 Москва, ул. Миклухо-Маклая, д. 6

А. А. Тютюнник

Российский университет дружбы народов

Email: tyutyunnik_aa@pfur.ru
Россия, 117198 Москва, ул. Миклухо-Маклая, д. 6

Д. А. Стариков

Российский университет дружбы народов

Email: starikov_da@pfur.ru
Россия, 117198 Москва, ул. Миклухо-Маклая, д. 6

Список литературы

Дополнительные файлы