Моделирование распространения поляризованного света в тонкоплёночной волноводной линзе

Полный текст

Введение В работе рассматривается задача математического моделирования дифракции электромагнитного поляризованного монохроматического излучения на трёхмерном утолщении волноводного слоя, формирующем структуру волноводной линзы на регулярном планарном диэлектрическом волноводе. Задачу можно рассматривать в декартовой системе координат, связанной с геометрией планарных волноводов. Плоские поверхности раздела между волноводным слоем, подложкой и покровным слоем параллельны (компланарны) плоскости xz. Ось y перпендикулярна этим плоским поверхностям. Излучение распространяется вдоль оси z (в направлении возрастания переменной z), геометрия системы до помещения в неё утолщения и невозмущённое электромагнитное поле излучения инвариантны относительно движений вдоль оси x, то есть ∂ ≡ 0. Тогда невозмущённое электромагнитное монохроматическое, поляризованное излучение удовлетворяет уравнениям Максвелла, которые в Гауссовой системе единиц в описанной системе координат распадается на две независимые подсистемы для TE- и TM-поляризаций. Если же нерегулярный участок интегрально-оптического волновода не обеспечивает инвариантности электромагнитного поля вдоль оси y, то на нем происходит деполяризация мод, т.е. все компоненты поля становятся отличными от тождественного нуля, что характерно для гибридных мод. Если нерегулярность мала, то и гибридизация (деполяризация) мала по сравнению с полем в регулярной части волновода [1-7]. Итак, рассмотрим распространение TE-мод в планарном регулярном волноводе, состоящем из трёх слоев: подложки, волноводного слоя и покровного слоя, характеризующихся различными коэффициентами преломления ns, n и n соответственно. Примем направление распространения волны за ось z, а ось y направим перпендикулярно слоям. Тогда Поместим сверху на основном волноводном слое небольшое утолщение - дополнительный волноводный слой переменной толщины h(x, z) порядка длины волны излучения на участке радиуса R порядка нескольких длин волн излучения (то есть за единицу измерения длины выбрана длина волны электромагнитного излучения), формирующий тонкоплёночную волноводную линзу [8]. Дополнительный волноводный слой имеет тот же показатель преломления, что и основной слой, сечение получившейся волноведущей системы изображено на рис. 1. Рис. 1. Плоский волновод, вставленный в ящик R × Ry Будем считать, что дополнительный волноводный слой представляет собой тело Ниже для примера используется где d - параметр, характеризующий «малость» утолщения. Наша цель - исследовать дифракцию света на этой линзе. Описание приближенной математической модели Отметим, что рассматриваемое локальное утолщение волноводного слоя не обеспечивает инвариантности поля вдоль оси y. Однако, в силу малости возмущения регулярности, эффекты деполяризации малы по сравнению с основным вкладом поля регулярного волновода. Это позволяет на начальном этапе рассмотрения пренебречь векторным характером распространяющегося электромагнитного излучения. В скалярном приближении распространение волны можно описать уравнением Гельмгольца а малая добавка dq1 характеризует утолщение. Поставить сразу парциальные условия изучения невозможно, поскольку спектральная задача для оператора ∆ + k2q0(x, y) на R2 имеет смешанный спектр: дискретный и непрерывный. Простейший путь к постановке условий - рассмотреть открытый волновод в закрытом волноводе, то есть ограничить рассматриваемое пространство R3 компактной областью G : x < R, y < Ry , изображённой на рис. 1 пунктиром. Мы полагаем, что в реальной системе объекты, помещённые достаточно далеко от волноводного слоя, не влияют существенным образом на интересующие нас характеристики волны. Данное предположение вносит дополнительное возмущение в задачу. Мы предполагаем его малым в начале рассмотрения. Результаты численных экспериментов, проведённых нами в последующих разделах, подтверждают сделанное предположение, поэтому примем u y= R = 0, а также, что поле на некотором удалении от утолщения не зависит от x, поэтому Остаётся поставить парциальные условия излучения. Обозначим собственные значения и собственные функции задачи как ln и n соответственно. Система собственных функций этой задачи полна в пространстве L2(G) [9]. Волну, падающую на линзу, при z < -L можно описать как разложение по собственным функциям, отвечающим отрицательным собственным значениям, причём ign = ln, N - число отрицательных собственных значений задачи (1). Отражённую от линзы волну при z < -L можно представить как а прошедшую при z > L как где Rn, n - неизвестные числа, а суммы распространяются на все собственные функции задачи (1). Задача отыскания функции u и последовательности чисел n , Rn , именуемых коэффициентами прохождения и отражения вида имеет (и притом единственное) решение и является корректной задачей математической физики [10, 11]. Мы будем использовать её в качестве математической модели для описания поля открытого волновода. Решение задачи дифракции в первом порядке теории возмущений Чтобы получить простые формулы для решения задачи (2), рассмотрим заполнение вида и будем искать решение в виде в рамках теории возмущений, возможность применения которой к задачам с парциальными условиями излучения была обоснована в [11, 12]. Для отыскания возмущённой части решения вида имеем соотношение вида или, после проектирования на n, вида 2 ∞ где Функция Грина для уравнения (3) с парциальными условиями излучения может быть выписана явно [13], а само решение с её помощью записано в виде Таким образом, полное поле даётся формулой где поправки первого порядка малости по d имеют вид: Следует заметить, что q1 имеет компактный носитель, поэтому на самом деле здесь интеграл распространяется на конечный отрезок. При z > L модуль |z - | = z - , поэтому Полагая можем переписать предыдущую формулу как Теорема 1. В первом порядке теории возмущений коэффициенты прохождения и отражения n-ой моды даются формулами где Матрицы d и dR будем называть матрицами прохождения и отражения. Границы применимости модели С точки зрения скалярной волновой оптики поле в открытом волноводе, помещённом в закрытый волновод, представляет собой линейную комбинацию волн где - собственные функ√ции задачи (1), отвечающие отрицательным собственным значениям, причём ign = ln. Собственные функции задачи (1) уже при небольших частотах k чётко распадаются на два класса: локализованные и не локализованные в волноводном слое. Первые мы будем интерпретировать как канализированные среднем слоем, а вторые - как канализированные объемлющим закрытым волноводом - ящиком. Гипотеза, лежащая в основе предложенной модели, может быть сформулирована следующим образом: если падающая волна представляет собой суперпозицию локализованных мод, то коэффициенты прохождения и отражения n и Rn локализованных мод не зависят заметным образом от параметров R, Ry ящика. Расхождения, которые неизбежно возникают при применении ящиков различных размеров, несут информацию, полезную для численного анализа модели: • n-ю моду можно вычислить при нескольких различных значениях R, Ry , величина dln характеризует дисперсию, которая не может быть уточнена в рамках рассматриваемой модели, вклады нелокализованных мод в любую величину характеризуют в своей совокупности рассеянную энергию, сами же по себе не имеют физического смысла. Так называемые «точные» модели также привносят погрешности, оценка которых зачастую не производится. Чем выше номер локализованной моды, тем более заметна зависимость собственных значений от R, Ry . Однако идея предельного перехода R плоха не только с вычислительной точки зрения. Это обстоятельство указывает на то, что рассматриваемая модель - модель с парциальным распределением точности [14]. Вычисление локализованных собственных функций Все собственные значения задачи (1) можно найти по методу разделения переменных. Полагая и подставляя (4) в (1), получим одномерную задачу на собственные значения Поэтому достаточно найти все собственные значения при s = 0, остальные получаются из них сдвигом на ((s)/(2R))2. При этом собственное значение l отвечает локализованной моде, если (y) экспоненциально убывает в покровном слое и подложке [15-17]: -k2n2 < l - ((s)/(2R)) < min(-k2n2, -k2n2). Отсюда, будет ли мода локализована при s > 0 или нет, существенно зависит от выбора R. Обратимся к вычислению локализованных мод, не зависящих от x (s = 0), сшиванием решений на двух разрывах коэффициента преломления. Для удобства опишем обозначения, которые используются в пакете Lu¨neburg под Sage. • 0 - подложка, 1 - волноводный слой, 2 - покровный слой, • n=[ns,nf,nc] - список с показателями преломления, • h=[h0,h1] - отрезок оси y, который занимает волноводный слой, Функция luneburg_eigenplot(n,h,k,Ry) проводит численный расчёт дисперсионной зависимости и строит график дисперсионной кривой (в логарифмическом масштабе), нули которой суть искомые собственные значения. Будем далее рассматривать размеры волноводной структуры в единицах длины волны электромагнитного излучения. Рассмотрим открытый волновод с толщиной волноводного слоя в 1 длину волны и помещённого в ящик с Ry = 10 длин волн, при k = 3 этот график можно построить командой: sage: load(’sage/luneburg.sage’) sage: var(’x,y,z’) sage: luneburg_eigenplot([1.1,2,1],[0, 1],3,10) На рис. 2 представлены графики при Ry = 2 и Ry = 10. Хорошо видно, что положение нулей не зависит от Ry . По графику видно, что имеется ровно две собственных функции, возле l = -31 и -18. Рис. 2. График левой части характеристического уравнения при n = [1.1, 2, 1], k = 3 для двух значений: (1) для Ry = 10; (2) для Ry = 2 Замечание 1. По графику хорошо видно, что задача имеет ту же вычислительную особенность, которая в [18] описывается в разделе, посвящённом построению графиков. Отделение корней здесь будет нетривиальной задачей. Для более точного вычисления корней и собственных функций служит функция luneburg_eigenfunction(n,h,k,Ry,lambda), где l - приближенное значение для собственного значения, найденное по графику, искомое собственное значение должно лежать на отрезке [l, l + 1]. Эта функция в качестве выходных данных предоставляет список, 0-м элемент которого служит уточнённое собственное значение, а 1-м элементом - собственная функция, описанная как кусочно-аналитическое - выражение и нормированная на L2( Ry , Ry ). Для рассматриваемого примера имеем два собственных значения: sage: luneburg_eigenfunction([1.1,2,1],[0, 1],3,10,-32) [-31.00587372275704, piecewise(y|-->-(1.749050320801987e-20)* e^(4.48507232079451*y + 44.8507232079451) + (1.749050320801987e-20)*e^(-4.48507232079451*y - 44.8507232079451) on (-10, 0), y|-->-0.5262885275712399* cos(2.234754187207837*y) - 1.0562423917910484* sin(2.234754187207837*y) on [0, 1], y|--> (2.3432737836960847e-19)*e^(4.69104185898581*y - 46.9104185898581) - (2.3432737836960847e-19)* e^(-4.69104185898581*y + 46.9104185898581) on (1, 10); y)] sage: luneburg_eigenfunction([1.1,2,1],[0, 1],3,10,-18) [-17.624762558080413, piecewise(y|-->(4.939108333233261e-12)* e^(2.59514210749246*y + 25.9514210749246) - (4.939108333233261e-12)*e^(-2.59514210749246*y - 25.9514210749246) on (-10, 0), y|-->0.920889502548748* cos(4.286634745569021*y) + 0.5575093905265337* sin(4.286634745569021*y) on [0, 1], y|--> (2.9481241586952706e-12)*e^(2.936794606042515*y - 29.36794606042515) - (2.9481241586952706e-12)* e^(-2.936794606042515*y + 29.36794606042515) on (1, 10); y)] Графики собственных функции построены на рис. 3 и 4. Следует обратить внимание на то, что, например, в покровном слое собственная функция имеет экспоненциально растущий член с амплитудой 1.7 10-20, который на отрезке 1 < y < 10 не оказывает заметного влияния на поведение собственной функции. Однако едва ли разумно отбрасывать этот член, ведь второе слагаемое имеет тот же порядок. -4 -2 2 4 -0.2 -0.4 -0.6 -0.8 -1 -1.2 Рис. 3. График собственной функции моды E0 при n = [1.1, 2, 1], k = 3 1 0.5 -4 -2 2 4 -0.5 -1 Рис. 4. График собственной функции моды E1 при n = [1.1, 2, 1], k = 3 Волноводная линза Обратимся теперь к дифракции на утолщении волноводного слоя. Пусть для примера волноводный слой имеет утолщение в форме полусферы, заполненной веществом с показателем преломления nl . В полярной системе координат Если n и m не зависят от x, то nm имеет вещественное значение, поскольку При этом для вычисления удобно пользоваться выражением В частности, Для вычисления этого интеграла можно воспользоваться стандартной подпрограммой Sage: sage: h1=1 sage: R=1 sage: delta=0.1 sage: nf=2 sage: nc=1 sage: ns=1.1 sage: nl=1.8 sage: v=luneburg_eigenfunction([ns,nf,nc],[0, 1],3,10,-32)[1] sage: numerical_integral(pi*(nl^2-nc^2)*(v.expression_at(2))^2* (R^2-(y-h1)^2/delta^2),h1,h1+delta) (0.0871817000736067, 9.679113074768287e-16) Вторая величина из числа выходных параметров numerical_integral, - ошибка вычисления интеграла. В этом примере стенка ящика была помещена на расстоянии десятка толщин волноводного слоя от этого слоя (Ry = 10), если же придвинуть стенку ящика на расстояние 1 толщины волноводного слоя (Ry = 2), то величина d11 почти не изменится: sage: v=luneburg_eigenfunction([ns,nf,nc],[0, 1],3,2,-32)[1] sage: numerical_integral(pi*(nl^2-nc^2)*(v.expression_at(1.5))^2* (R^2-(y-h1)^2/delta^2),h1,h1+delta) (0.08716079779015289, 9.67679245513252e-16) Именно это обстоятельство свидетельствует о применимости используемой модели: перемещение стенок ящика не влияет сколько-нибудь заметным образом на величину 11. Остальные элементы матриц d и dR для первых двух локализованных мод можно вычислить тем же путём. Так при d = 0.1 получается 0.0871817000736067 0.0367243736664542 0.0367243736664542 0.299853381650346 0.0181280522835082 0.0345146071574481 0.0345146071574481 0.0656532582915378 при d = 0.01 получается 0.000482453707056942 0.0361978913181138 0.000273495110082350 0.000482170019481738 0.000482170019481738 0.000850060161436371 d = (︂ 0.0116699996604377 0.000482453707056942 dR = (︂ при d = 0.001 получается 0.00120423307443047 4.96623379696062 × 10-6 4.96623379696062 × 10-6 0.00369188780662724 2.83607252631433 × 10- 4.96620475470580 × 10- 4.96620475470580 × 10-6 8.69624769266800 × 10-6 Можно видеть, что с уменьшением d матрицы коэффициентов dR и d стремятся к нулю. Заключение В работах по исследованию оптических волноводов обычно без дополнительного обоснования принимают, что поле внутри волноводного слоя не зависит от электромагнитных явлений, происходящих на расстоянии в несколько десятков длин волн от него. Это, разумеется, не вполне верно: некоторые характеристики поля внутри волновода зависят от этих явлений, но исследователей интересуют те параметры волноводного излучения, которые очень слабо зависят от окружения волновода. Данное рассуждение является основанием для формулировки математической модели (2) волноводного распространения поляризованного света в плавно-нерегулярном интегрально-оптическом волноводе. Предложенная модель не является точной, так как часть излучаемой открытым волноводом световой энергии канализируется бесконечно высокими стенками потенциала на границе объемлющего закрытого волновода. Достоинством модели является корректность формулируемой задачи. Для исследования её решений, их зависимости от граничных условий и от коэффициентов уравнения можно использовать широкий ассортимент строгих математических методов и апробированных компьютерных программ. Этим преимуществом можно воспользоваться для исследования того возмущения, которое вносится при помещении открытого волновода внутрь объемлющего закрытого волновода. Проведение численных экспериментов позволило показать, что с уменьшением дельта матрица коэффициентов отражения dR стремится к нулю, а матрица коэффициентов прохождения стремится к единичной матрице, то есть d 0. Причём обменные вклады, которым соответствуют внедиагональные элементы матриц, стремятся к нулю на порядок быстрее, чем диагональные члены. Сравнение полученных в данной работе результатов с результатами работы [19], в которой неполным методом Галёркина была исследована аналогичная волноводная конфигурация в случае закрытого волновода, то есть случай такого Ry , который приближает стенки объемлющего волновода вплотную к границам открытого волновода, показало качественное совпадение численных результатов.

Об авторах

Дмитрий Валентинович Диваков

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

Email: dmitriy.divakov@gmail.com
Кафедра прикладной информатики и теории вероятностей ул. Миклухо-Маклая, д. 6, Москва, Россия, 117198

Михаил Дмитриевич Малых

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

Email: malykhmd@yandex.ru
Кафедра прикладной информатики и теории вероятностей ул. Миклухо-Маклая, д. 6, Москва, Россия, 117198

Антон Леонидович Севастьянов

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

Email: alsevastyanov@gmail.com
Кафедра прикладной информатики и теории вероятностей ул. Миклухо-Маклая, д. 6, Москва, Россия, 117198

Леонид Антонович Севастьянов

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

Email: sevast@sci.pfu.edu.ru
Кафедра прикладной информатики и теории вероятностей ул. Миклухо-Маклая, д. 6, Москва, Россия, 117198

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