Модифицированный метод распространяющегося пучка и его применение к расчету распространения в волноводах с изменяющимся профилем показателя преломления Текст научной статьи по специальности «Математика»

МОДИФИЦИРОВАННЫЙ МЕТОД РАСПРОСТРАНЯЮЩЕГОСЯ ПУЧКА И ЕГО ПРИМЕНЕНИЕК РАСЧЕТУ РАСПРОСТРАНЕНИЯ В ВОЛНОВОДАХ С ИЗМЕНЯЮЩИМСЯ ПРОФИЛЕМ ПОКАЗАТЕЛЯ ПРЕЛОМЛЕНИЯ

А.В. Гаврилов

Самарский государственный аэрокосмический университет им. С.П. Королева Институт систем обработки изображений РАН

Аннотация

В работе рассматривается задача моделирования распространения электромагнитного излучения в волноводах с продольно-неоднородным распределением показателя преломления. В качестве базового метода моделирования предложен и рассмотрен метод распространяющегося пучка (Beam Propagation Method, BPM), выявлены его недостатки. Приведен вывод обобщенных однонаправленных уравнений Гельмгольца, позволяющих нивелировать ряд недостатков традиционного метода. Сформулирован модифицированный метод распространяющегося пучка, основанный на предложенных уравнениях, описана область его применимости.

Введение

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

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

Существуют различные методы, позволяющие моделировать распространение света в оптических элементах, однако не все они применимы к моделированию оптического волокна с изменяющимся поперечным распределением показателя преломления. Данная работа посвящена вопросам применимости метода распространяющегося пучка (Beam Propagation Method, BPM), широко применяющегося в волоконной и интегральной оптике, к моделированию процессов в оптических волноводах и модовому анализу распространяющихся пучков.

1. Анализ физических процессов в волноводах

с изменяющимся поперечным распределением показателя преломления

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

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

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

Особую роль в исследовании оптоволокна играет понятие моды - собственной функции оператора распространения, определяемого профилем показателя преломления [6]. Введенное в оптоволокно излучение в процессе распространения неизбежно (и обычно с потерей энергии) приходит к виду, когда пучок представляет собой моду или комбинацию мод, характерных для данного волокна [7]. Далее такие пучки распространяются без потерь энергии, а распределение энергии в тангенциальной плоскости определяется набором мод в пучке.

Основным средством внесения нелинейности в волновод является изменение распределения показателя преломления (часто достаточно существенное). Изменить профиль можно различными способами: изменением геометрии центральной линии волокна (например, микроизгибы [8]), изменением оптических свойств внешним воздействием (например, под действием давления [9]) или прямым изменением профиля (например, нанесением микрорельефов на торцы волокна [3]).

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

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

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

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

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

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

2. Метод распространяющегося пучка для волоконной и интегральной оптики

Изначально метод распространяющегося пучка был сформулирован Фейтом и Флекком в виде, отличном от современных модификаций, и основывался на послойном расчете распространения пучка излучения с помощью прямого и обратного преобразования Фурье [10, 11]. Основным достоинством этого метода является возможность расчета распространения на большие расстояния при относительно малых затратах памяти за счет того, что для вычисления характеристик электромагнитного поля на следующем слое необходимо знать их лишь на предыдущем слое.

В начале 90-х был предложен метод, получивший название конечно-разностного метода распространяющегося пучка (Finite-Differential Beam Propagation Method, FD-BPM), основанный на послойном конечно-разностном решении следствий уравнений Гельмгольца [12]. Математически дан-

ный метод отличается от предложенного Фейтом и Флекком, однако обладает тем же основным преимуществом: возможностью расчета распространения на большие расстояния. Благодаря этому данный метод и его разновидности [13] получил широкое распространение в волоконной и интегральной оптике. Тем не менее, метод обладает рядом особенностей и недостатков, для рассмотрения которых необходимо сначала понять физические основы данного метода.

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

V2 E + k02 n2 E = -V( E -V ln n2),

V2 H + kl n2 H = H

V2 n2

- niV (Vn2 - H),

(i)

где E и H - вектора комплексных амплитуд напряженности, соответственно, электрического и магнитного полей, к0 - волновое число распространяющейся волны в вакууме, n - показатель преломления среды. Такие уравнения являются следствием из системы уравнений Максвелла [4, 5] при следующих предположениях:

1) электромагнитные характеристики среды постоянны во времени,

jU (t) = const, s (t) = const, (2)

где s = s( r, t) и jU = u( r, t) - соответственно, относительные диэлектрическая и магнитная проницаемости;

2) отсутствуют источники поля,

P- 0, j - 0, (3)

где p = p(r,t) - плотность электрических зарядов, j = j (r, t) - вектор плотности тока;

3) среда является магнитно-однородной,

u( r ) = const. (4)

Нетрудно убедиться, что условия, описываемые выражениями (2)-(4), не накладывают ограничений, несовместимых с описанными ранее процессами в волноводах, поэтому уравнения вида (1) могут применяться для описания этих процессов.

Чаще всего FD-BPM формулируется для декартовой системы координат (хотя существуют и альтернативные подходы [14-16]), поэтому перепишем выражения (1), считая ось z основным направлением распространения, и представим их в матричной форме:

dz2

JL

dz2

E + AE = 0,

H + BH = 0,

где Е = ((, Еу, Ех) и Н = (Их, Иу, Иг ) - вектора

комплексных амплитуд компонентов электрического и магнитного полей, соответственно, а матрицы А и В являются матричными дифференциальными операторами и имеют следующий вид:

Лхх Л ^ ( В XX Ву В. ^

А = Лх Луу , В = Вх Ву Вуг . (6)

ч Лгх Лу Лгг ) ч Вгх Ву Вгг )

Компоненты этих операторов описывают взаимодействие компонентов полей:

д2 д2

Л„ =

Ау =

Л =

ЛУХ =

ЛУУ =

Л =

Лх =

йх2 д^+ дх 1

д 1п п2 д д21п п2

ду дх дхду

д 1п п2 д д21п п2

дг дх дхдг

д 1п п2 д д21п п2

дх ду дудх

д2 дх2 д2 ду2 + д 1п п2 ду ,

д 1п п2 д д21п п2

дг ду дудг

д 1п п2 д ---и. д21п п2

д 1пп2 д д21пп2 ,, 2

+-2— + к2 п

дх2 0

д21п п2

ду2

+ к02 п2

дх дг дгдх

д 1пп2 д д21пп2

Л =___|__

г ду дг дгду

д2 д2 д 1пп2 д д21пп2 ,2 2 =—т + —т +--+-2— + к02 п ,

дх ду дг дг дг

д2 д2 1 дп2 д ,2 2

Вх =—- + —- + —--+ к02п2 -

дх ду п дх дх

П гп2

Вху =-

В„ = ■

1

дп2 д д2п2 +

ду дх дхду

дп2 д д2 п2 +

Вх =-

1

дг дх дхдг

дп2 д д2п2 +

п у дх ду дудх

д2 д2 1 дп2 д ,2 2 '

В =—- + —- + —--+ к02п2 --

дх ду п ду ду

дп2 д д2п2 +

дг ду дудг

(7)

(8) (9)

(10) (11) (12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20) (21)

1 ( дп2 д д2 п2 В„ =-т| ——+

дх дг дгдх

1 (дп2 д д2п2 Ву =-т1 ——+

п I ду дг дгду

д2 д2 1 дп2 д ,2 2 V2 п2 В„ =—т + —- + ^--+ к02 п2--'¿г—

(22)

(23)

(24)

дх2 ду2 п2 дг дг

где оператор Vав означает дифференцирование только по координатам а и в .

Из выражений (13)-(15) и (22)-(24) следует, что для возможности послойного расчета поля необходимо исключить из рассмотрения продольные компоненты поля. Это можно сделать следующими способами.

1. Рассмотрение только ТЕ-поляризации, тогда

= 0, а выражение (5) для электрической составляющей примет вид:

д2

дг

2 Ег + АгЕг= 0,

(25)

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

2. Рассмотрение только ТМ-поляризации, тогда И = 0, а выражение (5) для электрической составляющей примет вид:

_д

дг

2 Нт+ ВгНг = 0,

(26)

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

3. Исключение из рассмотрения не компонент поля, но элементов матриц А и В , отвечающих за взаимодействие с продольными компонентами. Для этого, как видно из выражений (9), (12), (18) и (21), необходимо наложить некоторые дополнительные требования на среду распространения: потребовать «плавности изменения показателя преломления в направлении распространения» [12, 17]. При этом полученные из (6) уравнения в точности совпадут с (25) и (26).

Численное решение непосредственно уравнений (25) и (26) не дает искомого преимущества послойного расчета характеристик поля (поскольку либо требуется решение задачи во всем рассчитываемом объеме, либо решение не сходится), поэтому необходимо понизить порядок дифференцирования. Для этого преобразуем (25) и (26) следующим образом:

| - + ^ Нт= 0,

(27)

где г - комплексная мнимая единица, а корень из матричной функции определен на основе собственных чисел этой функции [18].

Скобки в выражении (27) определяют распространение волны в прямом и обратном направлениях. Если рассматривать только вперед-распространяющиеся волны, то получим однонаправленные уравнения Гельмгольца:

f+'^J E'=0, f+н —0,

(28)

откуда нетрудно получить математическую модель, используемую в FD-BPM:

d

dZ Ег=-'^,

dZ H—

(29)

На решении уравнений (29) и их следствий (в скалярных и полу-векторных случаях [13]) и основываются современные методы FD-BPM и FE-BPM (Finite-Element BPM). Сделанные в ходе получения выражений (29) предположения о характеристиках поля и среды делают данные методы вполне применимыми к задачам интегральной и волоконной оптики.

Однако для задач, связанных с расчетом распространения в волноводах с изменяющимся профилем показателя преломления, данный метод обладает рядом недостатков:

1) не учитывается отражение энергии в случае резкого изменения показателя преломления;

2) рассматриваются либо поляризации, либо присутствует требование малости изменения показателя преломления в направлении распространения;

3) отсутствует взаимосвязь между компонентами магнитного и электрического полей.

И если первый недостаток является неизбежной платой за возможность расчета распространения в существенно протяженных средах (впрочем, существуют итерационные модификации BPM, позволяющие частично решить и эту проблему [13, 19]), то вторые два не согласуются с условиями распространения в рассматриваемом типе волноводов. Например, требование «плавного изменения» показателя преломления означает близость к нулю элементов матрицы Axz и Ayz (что фактически означает

d ln nA ~ 0), а также B и B (что, в свою оче-

/oz xz yz

редь, означает 0n/gZ ~ 0). Нетрудно заметить, что

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

бенностями моделирования рассматриваемых волноводов.

Модифицируем метод таким образом, чтобы указанные недостатки были преодолены. Для этого необходимо изменить саму используемую модель (28).

3. Обобщенные однонаправленные уравнения Гельмгольца

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

(Dz ! ^

E +

(Ar | АЛ

в [в v I )

E = 0,

где

E = ( Ex, Ey, Ez ,0 f

Dz

A* —

dz

— 0

0 —

( A

Axz v Ayz

d

dz ) 0 ^ 0

(30)

(31)

(32)

(33)

а в - квадратная матрица, состоящая из 0. Из (30) несложно получить выражение, аналогичное выражениям (29):

(Dz ! в^

в \в v I )

E — -'.

(АД А* ^

I

E,

(34)

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

Сначала докажем вспомогательное утверждение. Выражение

d

—ET—±i4AEt

(35)

по сути эквивалентно выражению (25), которое можно представить в эквивалентной форме:

(36)

(Dz ! 2 'Ат ! в^

в Гв E + в Гв E — 0,

(7 (7 V 1 ) (7 | (7 V 1 )

что, в свою очередь, эквивалентно

E — ±i.

(Di_ 9

I

(37)

Таким образом, выражения (35) и (37) тоже эквиваленты, переход между ними можно осуществить, не смотря на наличие корня из матричного оператора.

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

А =

тогда

Г-А +1 I -А. А

—Г "I" У

V

Г Аг | А. А

V I У

(38)

(39)

Используя тот факт, что оператор А дифференцируем, и используя следствия из определения функции от матрицы [18] разложим правую часть (39) в ряд Маклорена:

- »

I - А = Ха а '

(40)

где а е К. - коэффициенты разложения, а А' - степень матрицы.

Найдем общий вид матрицы А', для этого также представим ее в виде блочной матрицы:

А' =

г 1 ь,!2 о 3 ьТ "ь,

(41)

где ' Ь1 - соответствующие блоки 2х2 для степени I. Теперь найдем общий вид этих блоков:

Г1 ь,., !2ь,.,А Г1 ь !2ь АГ-Ат +1! -а*А

А 1+1 =

Г1',

Найдем более общее выражение для 2 ь ,, для этого, используя понятие обобщенной обратной матрицы, перепишем выражение из (42) в следующем виде:

(2ь,+1 -2 ь,). А-1 =-(-Ат+1)

введем обозначения ^ =2 ь1 А-1, д = -Ат +1, получим:

^+1 -£,=-д.

Нетрудно убедиться, что

=-(1 -д)(1 -д)-1 =

1=0

= -(1 -(-Ат+1 ))(1 -1 + Ат)-1 = -(I-(I - Аг)) А-1,

тогда окончательно с учетом (42) получаем для (41):

Г/ . I ^ /, . чЛ А

А 1 =

(-Ат+1 -(I - Ат)).

I

Полученное выражение позволяет преобразовать ряд (40):

'(-Ат+1 -(I - АтУ)А-'А. ^

=х

а

в Т

Г(-Ат+1

—в—п

V I

I

(,

-х<

I-(I-Ат)'Ув! в;; л"! *-1 *

V I

! А-1А*

откуда имеем:

рь'+, = ь'.(-Ат+1),

2 ь,+, = - ь,. А. +2 ь,,

3 ь+, =3 ь .(-Ат+1),

4 ь,+, =-3 ь,. А. +4 ь,,

=х<

-х<

(Ат+1)

I вА

V

Гв|

I

'Т

V I

(I - Ат)

Л'

'=0

гв I а

V I У

Г в в А .в] а-"АГ

V I т У

что с учетом значений оснований последовательностей

1 ь =- Ат+1,

2 ь1 =- А*,

ь1 =в,

ь = I,

дает:

1 ь! =(-Ат+1)) ь, =-(-Ат+1) А. +2 ь'-1,

ь, = в, ь, = I.

(42)

=х<

(-Ат +1)

! вА

-I—

У

Ха'

( I » Л

'¡Хар - Ат)

=0

в V I

Г в в А ^! "А-1А",

V I т у

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

'=0

'=0

=0

=0

=0

=0

Vi - A =

(I - A в

v 1J

v I J

в 4T=l i

(в л/I-(I - Ат)Л

в V J _

( Ат вНв! л/А: I-- +--1-'--- в в\ в J V 1

с =

D

A A

n2k0 У s0 V dy dx

ATA.

в в л в a;1 a*

V I : J

A n

Км\м0 V dy dx

(48)

получим обобщенные однонаправленные уравнения Гельмгольца:

—Е: = -i

dz т

X

А/2Е: + A: /2CH

- 12с

(А:_! в' в !"в + ( | — 1/ > в \ А/2 А* . (43) i-H:=-i dz : H:+ В/^ DE:

1 V 1 J в i в V 1 J

Подставив (43) в (39), и результат в (34), с учетом (31)-(33), а также доказанного в начале утверждения эквивалентности, перейдем обратно к рассмотрению только тангенциальных составляющих поля, и получим окончательно:

—ET = -i

dz т

J\ET+ A

- y2 ( Axz 0 Y Ez

V Ayz

0

0

(44)

Упростим (44) и аналогично запишем уравнение для магнитной составляющей поля:

- х Г у

A Er = —i

dz т

^ HT = —i dz

at e: + a;

bT h + в

- X

v AzE j

(BH '

VByzHz J

(45)

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

Для уменьшения количества неизвестных в уравнениях (45) воспользуемся следствиями из системы уравнений Максвелла [5]:

E.z =

И = -

k •( rx HT),

(46)

описывающими связь продольных и поперечных компонентов поля (здесь к - орт аппликаты). Выражения (46) также можно записать и в матричной форме:

(47)

Окончательно, подставив (47) в (45) и введя обозначения

(49)

4. Модифицированный метод распространяющегося пучка Аналогично тому, как на однонаправленных уравнениях Гельмгольца (29) основывается BPM, опираясь на уравнения (49), можно сформулировать основные формулы для модифицированного BPM.

Т.к. с вычислительной точки зрения FD-BPM представляет собой конечно-разностное решение набора дифференциальных уравнений, разумно привести уравнения (49) к виду, пригодному для конечно-разностного решения. Основной проблемой в данном смысле является наличие степенных функций от матричных дифференциальных операторов: нахождение их значений по определению [13] весьма затруднительно, поэтому значения этих функций приближают другими функциями, более удобными в вычислительном смысле. Конкретный вид этого приближения определяет разновидность и некоторые характеристики метода.

Метод для малых углов распространения Как в начальной версии метода BPM [11], так и в ранних версиях FD-BPM применялось т.н. «приближение плавной огибающей» (Slowly varying envelope approximation, SVEA), позволяющее за счет введения некоторых дополнительных физических предположений понизить порядок дифференциального уравнения. С математической же точки зрения, SVEA-BPM представляет собой разновидность метода, получаемую аппроксимацией корня с помощью ряда Тейлора [13].

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

X * « 1 + 1 X, 2 2

Д -1X. 2 2

X " 12

Тогда (49) примут вид:

A

(50)

3 ^ (1 1

—E: = -i

dz : _V 2 2

д „ Г(1 1

—H: = -i I - + -

dz : V 2 2

Ет+f 2 - -2 А

т V 2 2 1

снт

вт I H +

1

Вт IDE:

Переход от матричной формы даст четыре уравнения (по одному для каждой из тангенциальных компонент) в частных производных, имеющих 1-й порядок по г . При этом по тангенциальным координатам первые два будут иметь 2-й порядок для компонент электрического поля и 4-й для компонент магнитного, а вторые два, наоборот, 2-й порядок для компонент магнитного поля и 4-й для компонент электрического. При этом, как следует из (9), (12), (18) и (21), правые части (51) не содержат производных комплексных амплитуд по г . Такой вид уравнений, описывающих распространение в среде, позволяет применить для их решения конечно-разностные методы.

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

Метод для больших углов распространения

В том случае, если основное направление распространение находится под существенным углом к предполагаемой оси распространения г , требуется применение более точной аппроксимации матричной функции. Чаще всего в этих целях применяются приближения Паде [20]. Применительно к случаю (49) они будут иметь вид:

±р? х

(1 - х)) = -

(1 - х )-X = 2-

X q1

j=0 p

(52)

(2) xj j

X q

j=0

(2) „j

где [p, q] - порядки аппроксимации Паде, а p(° е К и q(j) е К - коэффициенты аппроксимации,

определяемые из равенства построенной дроби разложению исходной функции в ряд Тейлора вплоть до члена p + q . С этой точки зрения, в частности, приближение (51) было получено при помощи аппроксимации Паде [1,0].

Такой подход позволяет приблизиться к исходной функции с точностью до большего порядка производной, что в данном случае и определяет максимально допустимое отклонение основного направления распространения от оси z .

Обычно [13, 21] в FD-BPM и FE-BPM используются аппроксимации порядков [2,2], [3,2] и [3,3].

Применение таких порядков к выражениям (49) возможно, однако (в силу наличия двух различных знаменателей после замены в правой части) приводит к достаточно высоким порядкам дифференциальных уравнений. Но при этом порядок по z также, как и для (51), остается равным 1. Таким образом, применение аппроксимаций Паде приведет к достаточно сложным шаблонам и СЛАУ при численном решении.

Замечания о методе

Наряду с FD-BPM, в последнее время широко применяется и конечно-элементная версия метода FE-BPM [17]: данный подход совмещает дискретизацию в направлении распространения и аппроксимацию оператора. К сожалению, форма уравнений (49) не позволяет очевидно построить выражения, аналогичные используемым в FE-BPM.

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

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

Заключение

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

Традиционно используемый в интегральной и волоконной оптике метод FD-BPM, основанный на решении однонаправленных уравнений Гельмголь-ца, не во всех случаях пригоден для моделирования распространения в рассматриваемых волноводах по причине дополнительно накладываемых им ограничений на характер профиля показателя преломления.

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

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

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

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

Работа выполнена при поддержке российско-американской программы «Фундаментальные исследования и высшее образование» (RUX0-014-SA-06) и РФФИ (гранты №07-02-12134, №06-08-01024, №07-07-00210).

Литература

1. Reinhard M. Integrated Optics, design and modeling //. Boston, London: Artech House, 1994.

2. Takahashi H. et al. Arrayed waveguide grating for wavelength division multi/demultiplexer with nanometre resolution // Electron. Lett., 1990. Vol. 26 (2). PP. 87-88.

3. Gavrilov A. V., Karpeev S.V., Kazanskiy N.L., Pave-lyev V.S., Duparre M., Luedge B., Schroeter S. Selective excitation of step-index fiber modes // SPIE Proceedings, 2007. Vol. 6605. PP. 660508-1-660508-6

4. Никольский В.В., Никольская Т.И. Электродинамика и распространение радиоволн: учеб. пособие для вузов // М.: Наука, 1989.

5. Ильинский А.С., Кравцов В.В., Свешников А.Г. Математические модели электродинамики: учеб. пособие для вузов // М.: Высш. шк., 1991.

6. Снайдер А., Лав Д. Теория оптических волноводов // М.: Радио и связь, 1987.

7. Хонина С.Н Многомодовые лазерные поля в ступенчатых оптических волокнах // Самара: Самар. гос. аэрокосм. ун-т, 2004.

8. Karpeev S.V., Pavelyev V.S., Khonina S.N., Kazanskiy N.L., Gavrilov A.V., Eropolov V.A. Fibre sensors

based on transverse mode selection // Journal of Modern Optics, 2007. Vol. 54 (6). PP. 833-844.

9. Gottlieb M., Chapter 2.3 // CRC Handbook of Laser Science and Technology, Vol. IV, Part 2 ed. M.J. Weber, Boca Raton: CRC press, 1986.

10. Feit M.D., Fleck J.A. Light Propagation in Graded-Index Optical Fibers // Appl. Opt., 1978. Vol. 17 (24). PP. 3990-3998.

11. Okoshi T., Kitazawa S. The Beam Propagation Method // Analysis methods for electromagnetic wave problems. Editor E. Yamashita, Artech House, 1990. Chapter 10.

12. Huang W. et al. The Finite-Difference Vector Beam Propagation Method: Analysis and Assesment // J. of Lightwave Technology, 1992. Vol. 10 (3).

13. Lu Y.Y. Some Techniques for Computing Wave Propagation in Optical Waveguides // Communications in Computational Physics, 2006. Vol. 1. PP. 1056-1075.

14. Heiblum M., Harris J. Analysis of curved optical waveguides by conformal transformation // IEEE Journal of Quantum Electronics, 1975. Vol. 11 (2). PP. 75- 83.

15. Plaum B., Wagner D., Kasparek W., Thumm W.M. Optimization of waveguide bends and bent mode converters using a genetic algorithm // Infrared and Millimeter Waves Conference Digest, 12-15 Sept. 2000.

16. Belousov A.A., Gavrilov A.V., Degtyaryov A.A. Numerical Solution of the Problem of Electromagnetic Radiation Propagation in Radially Symmetric Waveguide and the Error Analysis of the Solution // Optical Memory & Neural Networks, 2005. Vol.14 (3).

17. Chui S.L., Lu Y.Y. A wide-angle full-vector beam propagation method based on ADI preconditioner // Journal of the Optical Society of America A, 2004. Vol. 21. PP. 420-425.

18. Голуб Дж., Ван Лоун Ч. Матричные вычисления // М.:Мир, 1999.

19. Rao H., Scarmozzino R., Osgood R.M. A bidirectional beam propagation method for multiple dielectric interfaces // IEEE Photonics Technology Letters, 1999. Vol. 11 (7). PP. 830-832.

20. Baker G.A., Graves-Morris P. Pade Approximant // Cambridge U.P., 1996.

21. Stoffer R., Bollerman P.A.A.J., Hoekstra H.J.W.M., van Groessen E.W.C., van Beckum F.P.H. New true fourth-order accurate scalar beam propagation methods for both TE ans TM polarizations // Optical and Quantum Electronics, 2000. Vol. 32.