Расчет векторных мод оптического волновода Текст научной статьи по специальности «Физика»

РАСЧЕТ ВЕКТОРНЫХ МОД ОПТИЧЕСКОГО ВОЛНОВОДА

Котляр В.В, Шуюпова Я. О. Институт систем обработки изображений РАН, Самарский государственный аэрокосмический университет

Аннотация

В работе рассмотрен метод расчета собственных электромагнитных мод волноводов, основанный на методе согласованных синусоидальных мод. Соответствующее программное обеспечение, разработанное в среде Ма^аЪ, применено для расчета мод слабонаправляющего волокна со ступенчатым профилем показателя преломления.

Введение

Существует множество различных методов расчета собственных мод волноводов, но далеко не все из них могут обеспечить полностью векторный анализ [1-6], многие ограничиваются лишь приближенными результатами в рамках скалярной теории. Тем не менее, в литературе не однократно подчеркивалась необходимость именно векторного подхода к исследованию, например, таких сложных, но перспективных структур как фотонные волноводы.

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

Достоинство выбранного метода по сравнению с методами, основанными на конечно-разностном подходе [5, 6], состоит в том, что он позволяет получать в качестве результата непрерывные функции, также, как и метод [4], использующий функции Грина. Но подход, рассмотренный в работе [4] в силу своих особенностей может быть применен лишь к круглым волокнам, тогда как метод согласованных синусоидальных мод [1-3] - к волноводам с произвольным сечением.

Данная работа является логическим продолжением работы [1], в которой рассматривался скалярный случай. Также как в [1] здесь справедливы предположения о возможности разбиения сечения волновода на конечное число строк и столбцов, таким образом, что ни одна ячейка этого разбиения не содержала бы неоднородностей. На границах сечения предполагается наличие так называемых электрических или магнитных «стенок», определяющих граничные условия задачи. Основываясь на этих допущениях, строится алгоритм численного расчета констант распространения мод и построения функций их электрических и магнитных компонент.

Отличия векторного описания мод оптического волновода от скалярного

Принципиальное отличие векторного случая от скалярного в рамках метода согласованных синусоидальных мод состоит в том, что необходимо рассматривать локальные моды двух различных поляризаций - ТЕ и ТМ, поскольку обе они вносят свой вклад в формирование гибридной моды волновода, выражение для которой принимает вид:

^ у) = 1+¿^(^(у)], (1)

р=е,кк=1

где ^ представляет собой любую из электрических или магнитных компонент поля моды, а внешняя сумма соответствует суммированию по поляризациям: ТЕ- р = к, ТМ- р = е . Обычно, ТМ-поляризация характеризуется отсутствием продольной составляющей электрического поля Ег = 0, а ТЕ-поляризация - отсутствуем магнитной составляющей Бг = 0 . В отечественной литературе - наоборот. Рассмотрим подробнее локальные моды, входящие в выражение (1). Локальная х - мода, имеет вид:

( , V

и % х) =

к%}

$рк

ССБ^ЧX - X>)] +

и<арк )

к О)

"-хрк

(2)

БШ!

[к«( X - х(т))],

здесь учитывается поляризация, и,

лярного случая, введен множитель

в отличие от ска-

( . \2

, необ-

к

о

к(т) кРк /

ходимость которого будет объяснена далее. Величины и , входящие в выражение (1), оп-

$р ар

ределяются через параметры волновода и локальные у - моды, как показано в таблице 1.

Таким образом, выражение (1) принимает вполне определенную форму для каждой из компонент векторов Е (х, у) и сБ (х, у):

(т) Ех (х, у):

у (т)() ^

Е икк (х) <

(т)

к=1

/ \

(т) „

Еу (х у) = - ^ ' к=1

(т)

ек

(х)

у

(у) - ^ к=1

(т)

,(т) ( ) иек (х)

.(т). )

2 (т)

(у)

(3а)

ек

\ о У

(т)

Чк (у)

(т)

( ( у )

Е(т) (х, у)=-!икт)(х>

к=1

сБт)(х, у)}(х)

ФФт)(у)

Ффт)(у)

к=1

(х)

к

ко2 ^(т)(у):

к=1

к=1

(х)-

-чкт)(у),

(3б)

(3в)

(3г)

2

у

0

к

о

cBym \ x, y) = £

( m ) 'hk

( k (m) A

(x)

k=1

hk

фт)(y), (3д) cB(m(x,y) =-£«hm)(x)kфт)(У)-b'-íx. (3е)

k k=i k0 k=i ko

-улу)

Таблица 1. Симметричные и антисимметричные у -компоненты поля, выраженные через локальные моды

2

k

o

ТЕ ТМ

F F(m) sh F(m) ah F(m) se F(m) ae

Ex (К1 V k0 y) 0 0 V (m)( y) k|?£(m)( y)

Ey 0 0 ' k(m) kek k0 I V(m)(y) J ^(m)(y)

Ez 0 ^(m)( y) k0 (-kz ^ k 2 V k0 y V (m)( y) s(m)( y) 0

cBx 0 Ф( y) k 2 k 0 Г *z A V k0 y V ( m ) ( y ) 0

cBy ( к (m) ^ khk k 0 V 0 y 2 Ф^Ч y) 0 0 0

Это гибридные моды, из которых нетрудно получить по отдельности ТЕ- и ТМ-моды.

В таблице 1 и выражениях (3) ф(у) - рассмотренные подробно в скалярном случае [1] локальные у - моды, в данном контексте соответствующие ТЕ-поляризации, что означает абсолютную идентичность выкладок для этих мод с выкладками, проведенными для скалярного случая. Моды у (у), представляющие ТМ-поляризацию, имеют следующий вид в столбце т разбиения:

уЛу) = У1) ^^(у - у^+к^т^Су - у(п))]. (4)

,(n,l)

,(.n,l)

Здесь величины у/^' и у^'"' - неизвестные, подлежащие определению значения функции и производной локальной у - моды соответственно на границах разделения однородных областей (ячеек). Для них справедливы соотношения:

V

(n,r )

V

(n+1,l)

,(n)

,(n+1)

(5)

(6)

е

Для мод у(у) матричные равенства, вытекающие из требования их непрерывности (5) и условия (6), налагаемого на их производные, имеют форму:

(n'.r) s

■(n'.r)

= p(rí)w(n)P(rí-1)W(n-1) • • •W(2)P1

V

V

(1,l) s (1,l)

(7)

vT1)

vn+1l)

= Qf>i+1)y(>i+1)^+2)^+2) ,, y(N-1)QN)

Vfr)

VNr

. (8)

Здесь Р(г) и Q(г) также, как в скалярном случае матрицы вида:

p(') =

cosí

[ky°d(0] sin[ky°d(-)]/ k(y

- sin[ky-)d (í)]ky° cos[k y-) d(-)]

eo =

008

[k<° d(0] - sin[ky-) d(0]/k?

sin[ky-) d(0]k® cos[ky-) d(

y J y y

('+1) _

W

y (') =

0

,(-+1)

- = 1,N -1.

0

,(-+1)

- = 1,N -1.

(9)

(10)

(11)

(12)

Очевидно, что матрицы W (г+ ) и V « являются взаимно обратными и зависят исключительно от диэлектрической структуры сечения волновода. Соотношения (7) и (8) также, как их аналоги для ТЕ-поляризации, используются для расчета величин, входящих в характеристическое уравнение относительно к

2 .

е(п,)у{п'г)ч(ап'+и) - е(п'+1 )у(а"'г)у{"'+и) = 0 . (13)

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

Пусть кк - набор решений уравнения соответствующий ТЕ-случаю, а к1к - решения уравнения (13), описывающие локальные ТМ-моды.

Также незначительные изменения для локальных ТМ мод необходимо внести в формулу для расчета интеграла нормировки:

1 ' 1- [

) ..(1) Л

у - ч (у)

у(" +1 - у(1) > ((у)

йу = 1.

(14)

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

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

- затем рассчитать константы ф,("'1), фО"'1), Ф(п,г), ф(0П,Г), исходя из условия нормировки, чтобы окончательно сформировать К у - мод, соответствующих случаю ТЕ-поляризации;

- для каждого столбца сечения необходимо найти ровно столько же К корней ке2к уравнения (13), определяющих моды ч(у) с граничными условиями противоположными ТЕ-случаю. То есть, если для мод фу) использовались граничные условия, соответствующие электрическим стенкам, то для Ч(у) нужно использовать магнитные стенки и наоборот;

- далее рассчитать константы Ч"'1),

Ч((П,Г), Ча ',г), исходя из условия нормировки (14), чтобы окончательно сформировать К у - мод вида (4), соответствующих случаю ТМ-поляризации.

Существенной модификации требует и алгоритм поиска констант распространения и собственно «сшивки» локальных мод в непрерывные функции, описывающие компоненты векторного поля. Для х - локальных мод обеих поляризаций справедливы выражения:

и

(т,1) _

ар

= -Т.

(т)и(т, I) р -р

(т)и(т,1)

+ £

(т)ТТ(т,г)

и(т,Г) = О ( т)ит(т,1) + Т

ар р ,р р

и

р зр •>

( т)и( т,г)

(15а) (15б)

Здесь Т(рт) и Орт) - диагональные матрицы, содержащие следующие элементы:

трт = (ко/к(рт))2к(т/гап(к^й(хт)), (16а) =(ко/к(Рт))2к(т/*т(к<тк)йхт)). (Ш)

Равенства, связывающие значения локальных х-мод обеих поляризаций из соседних столбцов и их производные, приведены ниже:

тт(т,г) _ ^ (тт+^т-Дт+М) и к = °кк и-

зк

(17а)

и(атг) = окг+^иат1^ -кокт^чг10, (17б)

и-т,г) = о (т,т+1) и (т+1,1)

и(

(17в)

и(т) = оС^ит+У) + кртХт)и{тХГ), (17г)

где матрицы О^, Ое(ет,т,) и О^ имеют следующие элементы:

о^т =<Фкт)\Ф(т]) >,

кккр

(т,т,) =< ;(т) \ Ч (т') / ((т')/

(18а) (18б)

(т)2

О'й;> =< Фк \Чр /^'(у) > /ккк + < Ф(т)\Ч рт,)/((т0(у) > / к(;,)2. (18в)

Выражение для < фк \ фр > через значения функций и производных на границах слоев (кон-

станты

4("Л ф ("Л ф (п,г) ф (П,Г)

ка ' гз > га

) уже определялось ранее, в рамках скалярного подхода [1]. Далее приведем формулы удобные для расчета элементов остальных матриц интегралов перекрытия:

<Чк \ Чр / ((у) >=

-ЧкЧ^-ЧЧ +ч0?<))

п=1

(п)(4 - кр)

< Фк \ Чр/((у) >=

£(фЧ -кУ^ + фПг)ЧГ) + кипг)4г))

(п)(к2кк -кр

< фк \Чр/((у) >=

1 (п)(к2кк -кр

(19)

(20)

(21)

Нужно отметить, что введение множителя

С \2

к0

к

(т)

отличающего выражение для х - мод в век-

V рк )

торном случае от скалярного аналога, позволяет, сохранить простое выражения для нормировки (14) локальных у - мод ТМ - поляризации и получать выражения для элементов матриц Ок^) и Ое^), которые не зависят от корней характеристических уравнений к2рк.

Аналогично задаче отыскания констант распространения для скалярного случая, можно сформулировать задачу решения однородной алгебраической системы уравнений:

1(кг )и = 0,

(22)

где матрица Н(к2) имеет ленточную структуру и зависит от параметра кг:

2(кг)=

А? С (2) к В2 О О О О О О О

С™ А((2) С в(2) О О О О О О

^ О Ак(3) с (3) к Вк(3) О О О О О

О сез) А((3) С В(3 О О О О

О О О О

О О О О О О п (М ) к О Ак с (М) к

О О О О О О О п (М ) е с(м) А(М )

(23)

Компоненты блочной матрицы ), выражаются через уже определенные матрицы интегралов перекрытия и диагональные матрицы Т(р

'р и ер •

Ак = Скк Т И Скк + Тк , д(т) = 0(т-1,т)Т т'(т-1) 0(т-1,т) +

в случае магнитных стенок выражения для А(м) отличаются от общей схемы:

(24а) (24б) А™ и

42) = оИ2л)(Ти(1) -ек])(Тк°})-1 ек1))0кк2) + (24в)

а(2) = о£.2)Т(Те(1) -е(1)(Те(Ц)-1 8<Р)С™ + т(2\ (24г)

АМ =сИМм-1)тИм-1)сИМАМ) +ТМ) -е^Т^ГеМ, (24д) =сСМ-ХМ)Тт<М-11)сСМ-1М) +ТМ) ^(Т^)-^, (24ж)

5

(т

(т

5

<

С (т Ск

С (т Се

0ст

= - ек скк

= - е(т) С

I (т,т+1)

ее -

(т,т-1) (т-1,т) = К2Скк Ске

— Ь Г)(т-1,т)Т (т,т-1)Т = КЕСее Ске

(т,т-1) (т-1) = -Скк ек ,

= -С

(т-1,т)Т е(т-1)

(24е) (24 з) (24к) (24л) (24м) (24н)

В выражениях (24) индекс т изменяется в диапазоне 2, М •

В свою очередь вектор и в (22), также содержит значения на границах локальных мод как ТЕ, так и ТМ - поляризации.

~и (2,1) " и як и(2,1) яе

иЦ,')

и ал

и =

и (М ,1) и як (М ,1)

(25)

и

Задача (22) решается тем же самым способом, который был предложен для скалярного случая в работе [1]. В результате, после получения константы распространения и доопределения локальных х - мод, из полученного набора функций по форму-

лам (3) рассчитываются компоненты векторного поля Ех, Еу, Ег, сВх, сВу и сВг.

Расчет векторным методом основной моды слабонаправляющего оптического волновода Алгоритм расчета векторных мод волноводов, основанный на методе согласованных синусоидальных мод, был реализован с помощью пакета МаНаЬ 6.0.

В качестве тестового примера для проверки работоспособности программной реализации использовалась модель слабонаправляющего круглого волокна с показателем преломления в сердечнике псо = 1,47 и пс1 = 1,463 - в оболочке. Структура сечения данной модели изображена на рис. 1.

Рис. 1. Модель поперечного сечения круглого слабонаправляющего волновода, темно-серым цветом показаны области со значением показателя преломления псо = 1,47 ,

светло-серым - с пс1

= 1,463

Для длины волны излучения X0 = 1,3 мкм был произведен расчет всех шести компонент электромагнитного поля основной моды, которая представляет наибольший интерес в силу того что, исследуемая модель являет собой образец одномодового волокна, поскольку условие многомодовости

V =

пСо + п1.

X 0

>> 1

(26)

не выполняется. Здесь р- характерный размер сердечника (радиус) равный в данном случае 3 мкм.

у = 1,472 +1,4632 • Ъмкм ^ 2 078 1,3 мкм

На рис. 2 показаны распределения интенсивности рассчитанных компонент.

а) Ех

б) Еу

в) Ег

г) сВх

д) сВу

е) сВ2

Рис. 2. Распределения интенсивности соответствующих электрических и магнитных компонент основной моды волокна к

с индексом = 1,4642 к0

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

Оценить масштаб каждой из составляющих помогает таблица интегралов от функции интенсивности по области сечения £ :

I(Р) = ЦР(х, у)|2 ёхёу.

Таблица 2. Значения интегралов от функции интенсивности компонент основной моды

Р I (Р)

Ех 0,3568

Еу 4,4577

Ег -0,0021

сВх 9,5642

сВу 0,7679

сВг -0,0102

Как следует из таблицы 2, наибольший вес имеет компонента сВх, тогда как обе продольные компоненты Е2 и сВ2 очень малы.

Полученные функции, описывающие поле основной моды, также подверглись проверке на удовлетворение уравнениям Максвелла:

д_

ду д_

дх

Е2 + = -к осВх

Ег + 1кгЕх = 1к0сВу

дт Еу +д" Ех = -'косВ2

дх ду

д 2 ~В2 + ¡к2Ву = 1к0££0И0с Ех

ду

д 2 ~В2 + ¡к2Вх =-'к0ее0М0с Еу

дх

■д Ву +д

дх ду

Ву + — Вх = ¡к 0ее0М0с 2 Е2

д д

—сВх + — сВу - гк.сВ. = 0,

дх х ду у 22

— Ех + — Еу - ¡к2Е2 = 0.

дх х ду у 22

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

' 0.0005

СКО{Ко(Е = -¡к0сВ) =

0.0017, 0.0064

СКО(ЯоГВ = ¡к 0ее0^0 с 2 Е) =

0.0004 0.0074 0.0053

(27)

СКО(ШуВ = 0) = 0.0094 , СКО(БпЕ = 0) = 0.0046 .

Следовательно, компоненты, полученные с помощью векторного метода согласованных синусоидальных мод, удовлетворяют уравнениям Максвелла

с точностью 0,01, то есть рассчитаны верно. Нарушение же радиальной симметрии компонент можно считать следствием использования граничных условий, налагаемых по прямоугольному контуру.

Заключение

В настоящей работе был рассмотрен векторный случай метода согласованных синусоидальных мод, его программная реализация применялась для расчета основной векторной моды круглого слабонаправляющего волокна. Проверка полученного результата с помощью подстановки его в уравнения Максвелла подтвердила правильность расчетов.

Благодарности

Работа выполнена в рамках российско-американской программы «Фундаментальные исследования и высшее образование» (BRHE), а также поддержана президентским грантом РФ НШ-1007.2003.01.

Литература

1. Котляр В.В., Шуюпова Я.О. Расчет пространственных мод оптических волноводов с неоднородным поперечным сечением методом согласованных синусоидальных мод // Компьютерная оптика, 2003. №. 25. P. 41-48.

2. Sudbo A.S. Film mode matching: a versatile method for mode film calculations in dielectric waveguides // Pure Appl. Opt., 1993. V. 2. P. 211-233.

3. Sudbo A.S. Improved formulation of the film mode matching method for mode film calculations in dielectric waveguides // Pure Appl. Opt., 1994. V. 3. P. 381-388.

4. Stagira S. Full vectorial analysis of cylindrical waveguides using Green functions. // Optics Communications, 2003. № 225. Р. 281-291.

5. Lusse P., Stuwe P., Schule J., Unger H. Analysis of vectorial mode fields in optical waveguides by a new finite difference method // J. Lighhtwave Techn., 1994. V. 12. N. 3. P. 487-493.

6. Hadley G.R., Smith R.E. Full-vector waveguide modeling using an iterative finite-difference method with transparent boundary conditions // J. Lighhtwave Techn., 1995. V. 13. N. 3. P. 465-469.

7. Rogge U., Pregla R. Method of lines for the analysis of dielectric waveguides // J. Lightwave Techn., 1993. V. 11. P. 2015-2020.

8. Sztefka G., Nogling H.P. Bidirectional eigenmode propagation for large refractive index steps // IEEE Photonics Techn. Left., 1993. V. 5. P. 554-557.

9. Rahman B.M.A., Davies J.B. Finite-Elements Solution of integrated optical waveguides // J. Lightwave Techn., 1984. V. 2. P. 682-687.

10. Koshiba M., Maruyama S., Hirayama K. A vector finite element method with the high-order mixed-interpolation-type triangular elements for optical waveguiding problems // J. Lightwave Techn., 1994. V. 12. N. 3. P. 495-502.

11. Lusse P., Stuwe P., Schule J., Unger H. Analysis of vectorial mode fields in optical waveguides by a new finite difference method // J. Lighhtwave Techn., 1994. V. 12. N. 3. P. 487-493.

12. Hadley G.R., Smith R.E. Full-vector waveguide modeling using an iterative finite-difference method with transparent boundary conditions // J. Lighhtwave Techn., 1995. V. 13. N. 3. P. 465-469.

13. Lin P.-L., Li B.-J. Semivectorial Helmholtz beam propagations by Lanczos reduction // IEEE J. Quant. Electr., 1993. V. 29. N. 8. P. 2385-2389.

14. Lee P.-C., Voges E. Three dimensional semi-vectorial wide-angle beam propagation method // J. Lighhtwave Techn., 1994. V. 12. N. 2. P. 215-224