Методы теории рассеяния для решения задач дифракционной оптики

МЕТОДЫ ТЕОРИИ РАССЕЯНИЯ ДЛЯ РЕШЕНИЯ ЗАДАЧ ДИФРАКЦИОННОЙ ОПТИКИ

С.Г. Волотовский, Н.Л. Казанский, С.И. Харитонов Институт систем обработки изображений РАН

Аннотация

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

Введение

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

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

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

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

1. между источником и дифракционным оптическим элементом,

2. подложки,

3. модуляции,

4. между областью модуляции и регистратором.

Необходимо найти значение векторов электрического и магнитного поля в области регистратора.

Рис. 1. Оптическая схема.

1. Формальная теория рассеяния для бивекторного электромагнитного поля

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

- = HW,

Е1"

Е2

Ез

Е4.

Уравнения Максвелла для бивекторного поля имеют вид:

г дW

- — = HW,

к д2

W =

где Н - матричный дифференциальный оператор, называемый также оператором Гамильтона, Е, Нг -компоненты электрического и магнитного поля вдоль координатных осей х1.

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

В координатном представлении оператор Гамильтона имеет следующий вид: " 0 Л]

Н =

В 0

Л =-

-1 д х1 0

к 2 0 дх

1 "д х1 0

к 2 0 дх

е- 0 0 е-1

-д х2 д х1 -д х2 д ,

-д 2 д 1 х х

-д х2 д ,

е 0 0 е

0 1 -1 0

0 -1 10

где х - декартовые координаты, е - диэлектрическая проницаемость среды.

Операторное уравнение (1) можно формально проинтегрировать:

w(/) = и (/,г0 )w(/0) =

= ехр (- г ((- /0 )н )) (/0),

где

<Ю (_ 1)" / ^-1 иI(/,/0)=1 + I^"-1 I^"-2 X

"=1 /"-1 /"-2

I^ I^0 [Я(„-!)Н("-2).Н( )Н()],

где /=к2.

2

2

X

//

1 10

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

Произвольную функцию из рассматриваемого линейного гильбертова пространства представим в виде линейной комбинации:

W = 1 I"

п

Набор функций будем называть волновой функцией бивекторного электромагнитного поля в I-представлении. Каждому абстрактному оператору в

данном базисе можно сопоставить матрицу Н"Ш

=Е НШ

т

Система уравнений Максвелла в Е-представлении записывается в виде:

' дТт

к д-=е нт (z) л (

к дz п

В общем случае система базисных функций не является счетной. В этом случае суммирование заменяется интегрированием.

Базисные векторы в координатном представлении имеют вид:

Xy1 y2m

m1

m2

где 8(г) - функция Дирака.

Бивекторное поле имеет вид:

W (г1, г 2 )=

+ш m=4 1 2 / \

= J I WyymXyi y2m (г1, x 2 VyVy 2

В дальнейшем для записи выражений будем использовать правило Эйнштейна, согласно которому по повторяющимся индексам производится суммирование или интегрирование.

Если в пространстве существуют два базиса, связанные соотношением

Ут = ^

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

л=

щ

где Н(^)Ч, Н(I )Ч - матричные элементы в О- и Е-

представлении.

Пусть бивекторное поле имеет следующий

вид:

в области 1

в области модуляции

W(г1, г2, z)= fn (z) Fn (г1, г2 ), в области 4

W(г1, г2, z)= vn (z) Vn (г1, г2 ).

Пусть связь между базисными векторами в F-, G- и ^-представлениях имеет вид

Vm =

Fm =

На границе области модуляции и области 4 волновая функция бивекторного поля в G-

представлении имеет вид vm (0). Запишем ее в F-

представлении:

fn

Поле в произвольной точке в области модуляции имеет вид:

fn

где эволюционная матрица удовлетворяет уравнению

i dum

-=nn

к дz

с начальными условиями ПШ (0) = 5 пт .

Бивекторное поле в Е-представлении на границе области модуляции и области 1 имеет вид:

fn (-a) = Um (-a) fm (0). То же самое в Q-представлении: qn (-a) = Sm (f ^ q) fm (-a).

Подставляя выражение в

Полученное выражение можно рассматривать как систему линейных алгебраических уравнений

относительно коэффициентов ч" (-а) и ук (0).

Выражение (22) можно переписать в эквивалентной форме

В данном случае интегрирование дифференциального уравнения для эволюционной матрицы проводилось в направлении противоположном направлению оси z. Рассмотрим метод, в котором интегрирование эволюционного уравнения проводится в направлении оси z.

Пусть ч" (-а) - волновая функция бивекторного поля в ^-представлении на границе области модуляции и области 1. Запишем ее в Е-представлении

I" (-а) = БШ (ч ^ I) Чт (-а). Поле в произвольной точке в области модуляции имеет вид

W

fn (z) = Un (z) fm (-a),

где эволюционная матрица удовлетворяет уравнению

8

m=1

i dunm

= НП

к дz

с начальными условиями ит (-а) = 5 пт .

Бивекторное поле в ^-представлении на границе области модуляции и области 1 имеет вид:

/п

То же самое в ^-представлении

Уп

Подставляя выражение

Полученное выражение можно рассматривать как систему линейных алгебраических уравнений

относительно коэффициентов чп (-а) и ук (0).

Выражение (27) можно переписать в следующем виде

БР

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

2. Ковариантная запись системы уравнений Максвелла в криволинейных координатах.

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

Для записи уравнения Максвелла в криволинейной системе координат введем тензоры

En, Dn, Нп, Bn

где

Dn = egnmEm , Bn = gnmHn

е - диэлектрическая проницаемость среды.

Система уравнений Максвелла в криволинейных координатах в ковариантной форме имеет вид

(rotH )n = -ikDn (rot E )n = -ikBn

где оператор (rotF)n представляется в виде

e

(rotF )n =

ijn ( dF

2y[g

dx1

dF±

dxJ

g =

det ),

gij - компоненты метрического тензора в криволинейной системе координат.

Рассмотрим криволинейную систему координат, в которой метрический тензор имеет вид

gl3 = g 23 = ^

g 33 = 1,

е д3Е} = -¡к (^']Н] + е д}Е3, (34) еу'д3И} =-к (^¡Е] + еу'д}И3,

2епт3

2епт3

Выражаем Е3, И3 из последних двух выражений, и, подставляя в первую пару уравнений, получаем:

s1] д з E y =-ik (Jg)g11H] -

-вд,

Snm3

2iky[g

s1] дзHj =-ik (jg)giiE]- -

-eyd.

Snm3 (дnEm -дmEn)' 2ikjg

Рассмотрим случай ортогональных криволинейных координат. В этом случае система уравнений Максвелла относительно четырех компонентных бивекторов имеет вид:

k

где

W =

Н =

д ZW = HW,

( Ei ^ E 2 Hi V H 2 J 0 A B 0

где

-1 (д ° yi 0

k2 V 0 д y 2

(- д y 2 д yi ^ (

V д y 2 дyi J V

1 (д yi 0

k2 V 0 д у 2

(- д y 2 д yi ^ -Г

V д y 2 дyi J V

0

0

Ш)

22 Г~\

g

0

22 Vg 0

/

& )-i

18) 0

/

E 0 0 e

-8

0 g 224s

22 4g 0

Л

/

Отметим, что при переходе от одной системы криволинейных координат к другой компоненты че-тырехкомпонентного бивектора - верхняя и нижняя

X

X

0

X

X

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

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

Базисные функции ^-пространственно-частотного представления имеют вид:

-^а1а2,1

с I -7 I 1 2

51; ехр ик 1а1 х +а 2 х

5ъ ехр (¡к 1а1 х1 +а2х2'

2

2

531 ехр (¡к 1а1 х1 +а2х 541 ехр (¡к 1а1 х1 +а2х

где 51к - символ Кронекера.

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

/ \ к=4 1 2, ,

Раа , х2 )=Е [ 5уук- ( - х)>

а1а2, г \ > / ¿-^ J а1а2, г ^ '

к=1

> Ху!У2к (/1' х2 к^ АУ 2.

>У, к1

Связь между волновыми функциями координатного и ^-пространственно-частотного представления имеет вид

к=4

WУ

к=1 -

= у Г 5 уа 2к.

J а1а2,г V >

где

х / а1а 2^а1 dа 2

(х1, х 2) 1 \ ' /

= о((1 - у1 )о(;

X 1 2 (х1, х2 )=

у у , т у

22

х 2 - 2

т1

т4.

Матрица преобразования имеет вид ^¿к

Матрица обратного преобразования имеет вид

^2к (- /)=

4п

ехр (-¡к (а1х1 +а 2 х2 ^О^

Гамильтониан в пространственно-частотном представлении имеет вид

Н°

12

В

0

^а

ла

1й 2

ОЪЮ

1Ш2

а1а2 ¡а1 0 "

л е-1

0 ¡а 2 _

- ¡ю2 ¡ю1 " 0 1"

X - -1 0

- ¡Ю 2 ¡ю1

1 (а1 -ю1;а2 - ю2)>

¡а1

0

0 ¡а 2

- т2 /ю1

- ¡ю /ю.

-е 1

0 1 -1 0

5Ю = 5(а - ю) - дельта-функция Дирака.

Уравнения Максвелла в пространственно-частотном представлении приобретают вид:

i д/а^а2г

= Н а^! / 12к

к дг

Для вакуума уравнения приобретают вид:

i д/(

= Н а^* I а1а2к к д2 а1а2к

ва

а1а 2 0

Л а1а2 =_

а1а2

¡а1 0

¡а 2

¡а 2

Г- ¡а 2 ¡а1 " 0 1"

|_- ¡а 2 ¡а1 -1 0

Г- ¡а 2 ¡а1 " 0 1"

|_- ¡а 2 ¡а1 -1 0

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

/ а

гахаг,\ 1 ±е /-01^2,2 1 ±е у-а1а2,3 1 ±е

= 1 Н

х ехрI ±

/±Ь =

/• аа,

±е

± ¡к ^ 1 -а12 -а^ ^

-а12 -а^ -а12 -а^

а.

а2

/•а1а2,1

1

1 ±е -а1а2,2

±е

•а^1а2,3

1

±е / а1

1 ±е

,-а1а2,4 1 ±е

=

-а2

а!

Л/1-а1

-а? -а2 1а!

-а12 - а2 |а2

х ехр I ± ¡к

: ¡к ^ 1 -а12

22 2 -а22

5а1 оа1

ЮЮ ю

1

В а1а2 =

5а1 оа1 -

юю

юю

2

11

0

Н а1а2 =

аа

12

12

0

0

0

В а1а2 =

аа

12

«

х

2

-а

2

X

12

0

г

1 +

22 1 -а, -а2

Введем К-пространственно-частотное представление (^-представление). В качестве базисных функций К-представления выберем собственные функции оператора Гамильтона для бивекторного поля, распространяющегося в вакууме. Базисные функции К-пространственно-частотного представления связаны с базисными функциями Е-представления

К

12п

2'

а.а 2п

= Б(V ^ /К^п (, X2)), сгп (V ^ / )=52 5»; Б"т (V ^ /),

Б"т (V ^ /) =

_г1п 1

/+е т=1

+е

/•1п г\

+н т = 2

1п е

/1п *

т = 4

/-е т = 3

Тензор обратного преобразования выглядит следующим образом

БЮ1Ю2т

а1а2т V / а2 а1 X /п т > V /

где матрица Р^ есть матрица эрмитово-сопряженная по отношению к матрице Бп (у ^ /),

а матрица

мп = Р^ (V ^ /). Здесь

верхний ин-

декс означает номер строки в матрице, а нижний индекс - номер столбца.

М =

1

0

ж ±

ж ± 0 0 ж ± 10

0 Ж ± 0

м - =

1

-ж ±

0

-ж ± 0

0 -Ж:

-ж ± 0 1

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

/ а1а2п = ба1а2п ( . /) ^ 12п V ^ /

а1а2п = Б а1а2п ( / . v) /12п 12п V./ У-'

4 Распространение бивекторного электромагнитного поля в вакууме

В настоящем разделе рассмотрим задачу распространения бивекторного электромагнитного поля в свободном пространстве.

Пусть в плоскости z=0 бивекторное поле имеет

вид:

ж = ж0 х1, х2 = ж

т хх к

(0).

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

/а1а2,(0) = Б*1*2-,! (х ^ /)жххк (0).

^^

Матричные элементы гамильтониана в вакуу-

ме в К-представлении И! дующей матрицы:

ю1ю2к

есть элементы сле-

Иа1<Х2 (V) = ^к 51 51 х

V > 1 Ю1

: - а^ - а

10 0 0 0 10 0 0 0 -10 0 0 0 -1

Решая систему уравнений Максвелла в К-представлении для произвольной плоскости, получим

V а1а21

V а1а22( z) v а1а2 3

Ч z)

V а1а2 4

V а1а21 (0) ехр ^ ¡ky|l—а[~-~а2z ^ Vа1а22 (0) ехр^ ¡кф -а2 - а2 z^

V а1а23

V а1а24(0)ехр ¡ку[[

- ¡кл11 - а,2 - а2 z

.

Далее перейдем от К-представления к координатному представлению

/ а1а2п (z) = ^^2« (V ^ /) V 12п (z) ,

жх

(z)=Б:а :к ( ^ х)/а1а 2 (

Выражение описывает бивекторное электромагнитное поле в плоскости регистратора.

5. Дифракция на пропускающих дифракционных оптических элементах

Рассмотрим дифракцию света на пропускающих дифракционных оптических элементах

Пусть Vа1<Х2п (-а) - волновая функция бивек-торного поля в К-пространственно-частотном представлении на границе области модуляции и области 1. Запишем ее в Е-пространственно-частотном представлении

/а1а2п

Поле в произвольной точке в области модуля-

ции имеет вид:

/ 12п (^ = и 12п (z) / О а1а2

(-а),

где эволюционная матрица ип (z) удовлетворяет интегро-дифференциальному уравнению

V

1

0

0

1

0

1

х

12

ди ю1ю2п

--а1°2т = Н

к д2 Р1Р21 у а1а2т >

с начальными условиями

ию

(-а) = 5Ю2 5Ю1 5п

Л ' г/., г,

Полученное выражение представляет собой интегро-дифференциальное уравнение.

Бивекторное поле в Б-представлении на границе области модуляции и области 4 имеет вид:

/ Юю2п (0) =

= иЮЮП(0)Б«^ ( - /)в1в2т (-а). То же самое в ^-представлении

V^ (0) = Б^т ( - V)/^ (0).

Подставляя выражение (72) в (73), получаем:

1а2П

в1 в2

Ю1Ю2 у

Б ^ (V - /) П2 к (-а).

Полученное выражение (74) представляет собой связь между полями V ^^ (-а) и Vа1 а2к (0) на границах области модуляции. Выражение (74) можно переписать в следующем виде:

Б ^ (V - /) а1а2: (0) =

= иЮЮР(0)БП1 !2к(V - /)п П2к(-а).

На практике поле удовлетворяет следующим условиям.

Условие 1. Поле прошедшее через оптический элемент не содержит волн, распространяющихся против оси 2:

V п П23(0) = V П24(0) = 0.

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

/а1а2,*

'(-а) = БЮЮ:( - V)/Ю1Ю2:(-а).

Решая систему интегральных уравнений (75) и используя результаты, полученные в предыдущем пункте, определяем прошедшее бивекторное поле в области 4.

6. Дискретное пространственно-частотное представление

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

пространственно-частотном представлениях имеет следующий вид:

/пт

, к

(0),

' d: ^ 2

d:d 2 - размеры дифракционного оптического элемента или периоды периодической двумерной структуры.

Тензор обратного преобразования имеет вид:

Бт,к

Гамильтониан в дискретном пространственно -частотном представлении имеет вид

Нпр =

т д

0

впр

т д

ЛпР т д

ЛпР =-д

/а1

'п-т, р-д

01 -1 0

о т о р

ЛпР = -

тд

- ¡а 2

¡а.

о т о р -

01 -1 0

'п-т, р-д •

Формулы (75) в дискретном пространственно -частотном представлении имеют вид:

Р (V - /(0) =

туд V

= итд (0)Б2д (V - /Уг (-а).

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

Бкдт (V - / )=о к орБт (V - /),

Будт (/ - v)=5ij ор (- )пк?кп.

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

I дитр

=Нпк

к д2 с начальными условиями

итур(-а)=о то)о р.

Формулы (84), (85), (86) могут быть получены из формул (56), (57), (70)< если заменить интегралы на интегральные суммы.

12

0

V

X

д

0

-/а

х

0

-/ а

0

- / а

0

12

V

7. Реализация вычислений Для того чтобы получить систему линейных алгебраических уравнений (83), необходимо многократное уравнение для эволюционной матрицы (86). Решение эволюционного уравнения сводится к решению системы обыкновенных дифференциальных уравнений с различными начальными условиями. Это позволяет использовать технику параллельных вычислений. В данном случае был использован под-

ход, состоящий в том, чтобы система уравнений для различных начальных условий решалась на различных компьютерах. Приводимые ниже результаты были получены с использованием кластера, состоящего из четырех двухпроцессорных компьютеров РЕЭТЩМ-П с частотой 450 МГц. Для решения системы обыкновенных дифференциальных уравнений использовались методы матричной экспоненты и Рунге-Кутта.

в

б

б

а

В

в

в

г

В

в

д

е

В

в

з

Рис. 2. Результаты расчета проекции вектора Умова-Пойтинга на оптическую ось для дифракции плоской электромагнитной волны на радиально-симметричной бинарной линзе в плоскости 2=0 (а, б), в плоскости 2=//2 (в, г), в плоскости 2=/(д, е), 2=3//2 (ж, з).

в

в

а б

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

В качестве примера был выбран расчет поля от бинарной цилиндрической линзы.

При расчете были использованы следующие параметры: радиус апертуры Л=4,82Х, фокусное расстояние /=4,82Х. На рисунках 2, 3 приведены результаты расчета проекции вектора Умова-Пойтинга на оптическую ось для дифракции плоской электромагнитной волны на бинарной ради-ально-симметричной линзе на различных расстояниях от оптического элемента. При этом на рисунках черный цвет соответствует максимальному значению. Приведенные выше результаты показывают работоспособность приведенного алгоритма. Разработанный метод снижает требования к вычислительным ресурсам, по сравнению с многомерными разностными методами. Анализ полученных результатов показывает, что максимум интенсивности в фокусе линзы примерно в четыре раза меньше значения, рассчитанного в скалярном приближении Кирхгофа. Однако, несмотря на это, бинарная линза сохраняет свои фокусирующие свойства. Следует отметить, что в отличие от скалярной теории, фокальное пятно имеет слегка вытянутую форму. Это объясняется тем, что в рамках электромагнитной теории не существует радиально-симметричных ре-

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

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

Литература

1. Ландау Л. Д., Лифшиц Е.М. Теоретическая физика // М., Наука. Т. 3. 1973. 504 с.

2. Методы компьютерной оптики // М., Физматлит. 2000. 688 с.