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

УДК 621.39

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

слоем

К. П. Ловецкий, А. А. Хохлов

Кафедра систем телекоммуникаций Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия

Моделирование взаимодействия ТЕ и ТМ поляризаций с оптическими структурами, характерные размеры которых могут быть сравнимы с длиной волны видимого диапазона, представляет большой практический интерес в связи с современными тенденциями в развитии микроэлектроники, оптоэлектроники и медицины. В США и Европе активно реализуются программы по развитию солнечной энергетики, основанной на использовании солнечных батарей, проектирование и оптимизация работы которых требует построения моделей точного учёта взаимодействия света с оптическими покрытиями. Эти модели также используются при создании и улучшении жидкокристаллических дисплеев и дисплеев с применением технологии OLED. В современных оптоэлектронных устройствах востребованы высококачественные зеркала, поляризаторы (общей толщиной < 1 мкм) и многие другие элементы дифракционной оптики. Также надо отметить использование многослойных структур в медицине. Из таких структур на основе кол-лагеновых плёнок изготавливают различные импланты, одним из примеров которых является искусственная роговица глаза.

Ключевые слова: анизотропная оптика, матрица Берремана, распространение электромагнитной волны, граничные условия электромагнитного поля, уравнения Максвелла, частотная дисперсия.

1. Введение

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

Будем рассматривать плоскую монохроматическую волну, характеризующуюся определёнными значениями частоты ш и волнового вектора к, удовлетворяющую однородному волновому уравнению, иначе называемую «нормальной» электромагнитной волной. Нормальные волны бывают разных видов. В рамках данной задачи будем рассматривать однородные электромагнитные волны, электрическое поле которых имеет вид [1]:

E = Eoe-i(wí-kf), k = -п(ш, s)s, (1)

где Ео — амплитуда поля — комплексный вектор, который не зависит от времени и координат г и ¿.В общем случае п = п + г\ — комплексный показатель преломления среды, а § = — единичный вещественный вектор. Анализ распространения

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

Статья поступила в редакцию 6 ноября 2009 г.

Работа выполнена при поддержке темплана РУДН, тема 020614-1-173.

Отметим, что в рамках данной задачи будет рассматриваться линейно поляризованная волна. В зависимости от направления её поляризации она всегда может быть рассмотрена как суперпозиция двух поляризаций — ТЕ моды и ТМ моды в выбранной декартовой системе координат.

2. Математическая модель

В качестве математической модели для решения этой задачи используем уравнения Максвелла в гауссовой системе единиц в отсутствие сторонних токов и зарядов:

rot Ё(г, i) = у B(r, i), rot Й(г, i) = - у D(f, i). (2)

Зависимость роторов полей от времени в случае волн вида (1) выражается следующим образом:

rot Ё(Г, i) = e-iwirot Ё(Г), rot ЙЙ(?, i) = e-iwirot Й(Г). (3)

Чтобы система уравнений была полной, её необходимо дополнить материальными уравнениями:

D(w,f) = е(ш)Ё(ш,г), В(ш,г) = Й(ш,г), (4)

где £ — тензор диэлектрической проницаемости, ш — циклическая частота, r — вектор, описывающий координаты полей.

Обоснуем принятый вид материальных уравнений. В рамках задачи рассматриваются немагнитные среды, это определяет тождественное равенство магнитной индукции и напряжённости магнитного поля. Теперь обоснуем зависимость диэлектрической проницаемости от частоты. Лишь в очень ограниченном круге задач можно считать компоненты тензора е постоянными. На самом деле оптические свойства материала зависят от частоты поля ш (эффект частотной дисперсии) и от волнового вектора к (эффект пространственной дисперсии), т.е. е = е(ш, к). Рассмотрим вкратце физическую природу частотной и пространственной дисперсии. Если собственные частоты среды ш8 меньше или сравнимы с частотой ш, то установление электрической поляризации не происходит мгновенно, т.е. поляризация поля, а значит и его индукция, определяется значениями поля не только в текущий момент времени, но ещё и в предыдущие моменты. Это и приводит к зависимости £ от ш при использовании Фурье-компонент полей [1]. Для большинства сред собственные частоты лежат в оптическом диапазоне, и поэтому частотная дисперсия в оптике должна быть учтена, т.е. е^ = е^(ш). Такое приближение является достаточно хорошим при решении практических задач, и оно будет использоваться в данной статье. Верхнее уравнение в системе (4) указывает на учёт частотной дисперсии среды. Однако это приближение является недостаточно точным. Например, в рамках него не могут быть объяснены некоторые известные явления, например гиротропия. Для более полного описания характеристик среды необходимо также учитывать зависимость тензора диэлектрической проницаемости от волнового вектора к: е ^ = е^-(ш, к). С физической точки зрения пространственная дисперсия означает зависимость электрической поляризации среды в точке от состояния электрического поля только в этой точке, хотя в реальности она определяется полем в некоторой окрестности этой точки. Это определяется зависимостью свойств среды от длины волны, или, что то же самое, от волнового вектора к, так как k = s, где Л — длина волны в среде [1]. Пространственная дисперсия определяется коэффициентом ^, где а — характерный размер среды (радиус «области влияния», молекулярного действия, и т.д.) Практический интерес в оптике представляют конденсированные среды, в частности жидкости и кристаллы, где а имеет порядок постоянной кристаллической решётки или размеров молекул, т.е. 10-9 — 10-10 м. Длина волны оптического диапазона имеет порядок 10-7 м, поэтому показатель пространственной дисперсии мал, и в задачах моделирования

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

Геометрия задачи отражена на рис. 1. Электромагнитная волна падает на непоглощающий анизотропный плоскопараллельный слой под некоторым углом «1 из изотропной непоглощающей среды с действительным коэффициентом преломления п\ и выходит в изотропный непоглощающий слой с действительным коэффициентом преломления «2. Рассмотрение непоглощающих сред обусловлено, в первую очередь, упрощением записи условий на границах раздела сред. В определённом классе реальных задач такое упрощение вполне приемлемо. Лабораторная система координат выбрана так, чтобы падение световой волны происходило в плоскости XOZ. Это обусловливает отсутствие у — компоненты волнового вектора падающей волны, т.е. кгупс = 0. Плоскопараллельный анизотропный слой расположен параллельно плоскости ХОУ. Считаем, что волна падает из верхней полуплоскости, отражается вверх в неё же и проходит в нижнюю полуплоскость. Исходя из этих соображений, ось 02 выбираем как направление, вдоль которого рассматривается эволюция процесса взаимодействия электромагнитной волны со слоем. Таким образом, при распространении волны вдоль выбранного направления будет меняться только компонента кг волнового вектора, а компонента кх будет оставаться неизменной.

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

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

0 Ь

Рис. 1. Падение электромагнитной волны на анизотропный слой

образом:

(еа 0 0\

е = А I 0 £ъ 0 | А-1, (5)

\ 0 0 £с)

где А — матрица поворота на углы Эйлера, которая представляет собой произведение матриц поворотов на углы прецессии, нутации и собственного вращения. В дальнейшем изложении будем полагать, что кристаллографическая система координат уже приведена к лабораторной, в которой и будут рассматриваться все величины и уравнения. Тензор диэлектрической проницаемости в лабораторной системе координат выглядит следующим образом:

11 £12 £13\

21 22 £ 23 | . (6)

£ 31 £32 £ 33/

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

После подстановки материальных уравнений (4) в уравнения Максвелла (2), учитывая вид волн (1) и раскрывая операции ротора, получим систему из шести уравнений: четырёх дифференциальных и двух алгебраических:

г с д ш дх

^^ дд ^^ ___ск%

— ТТП-х = £21Еш + £22Еу + £2зДг--

ш дг ш

--~ — Ну = £ЦЕх + £12Еу + £1зЕх

---Ну = £31 Ех + £ 32 Еу + £ззЕг

ш у у

(7)

г°д гс д скх скх

— ^ТЕу = Нx,--= Ну +--Ег , -

ш д ш д ш ш

Из третьего и шестого уравнений системы (7) можно выразить Ег и Нг через оставшиеся компоненты:

Ег = -

в 31 Еж + £32ЕУ + ^Н

зз

тт _С кх

Нг = -Еу

ш

(8)

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

д_

д

(Ех\

ЕУ

Нх \НУ)

ш

—А

/Ех\

ЕУ

Нх \НУ)

(9)

где матрица коэффициентов А выглядит следующим образом:

/

А

.£31 скх

£33 Ш 0

£32 скх £33 ш 0

-£ 21 +

\ £11

£23£31 £33 £13£31 £33

— £ 22 +

£23 £32 £33

+ ) ^

£12 —

£13£32 £33

-1

V Ш ) £33

1_\

0 §23 скх

0 _ £13

£23 скх £33 Ш £13 скх £33 ш

(10)

0

0

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

Решение системы дифференциальных уравнений с постоянными коэффициентами можно искать в виде матричной экспоненты (Н — толщина слоя, Ф =

(Е X ; Еу ; Нх; Ну ) ).

Ф(Н)=ехр(^ ^Н^)Ф(0) = Р (Н)Ф(0), (11)

где Р(Н) — так называемая матрица Берремана для однородной среды [2]. Начальными условиями системы является известная информация о падающей на слой электромагнитной волне.

Для нахождения матричной экспоненты предлагается воспользоваться преобразованием подобия:

Д = 5-1Бв, (12)

где 5-1 и 5 — матрицы с левыми и правыми собственными векторами матрицы Д, а Д — диагональная матрица, состоящая из собственных значений матрицы Д. Наиболее предпочтительным для численной реализации нахождения собственных значений и собственных векторов матрицы является метод вращений Якоби [3]. Для применения этого метода необходимо, чтобы матрица Д была диагонализуе-мой. Надо отметить, что не любая матрица может быть диагонализована, но всегда может быть диагонализована матрица с различными собственными значениями. В случае, если матрица Д не может быть диагонализована (например, в некоторых случаях при наличии кратных собственных значений), можно воспользоваться известной теоремой об аппроксимации [4]: для заданной матрицы А всегда можно найти матрицу В с различными собственными значениями, причём такую, что || А — В|| < 5, где 5 — любая наперёд заданная константа. Так как задача решается численно с применением компьютера, она решается с точностью до наперёд заданной машинной точности, и мы вправе применить эту теорему, используя машинную точность в качестве константы ё.

Известно, что в случае представления матрицы коэффициентов в виде (12) [5], функция от матрицы записывается следующим образом:

ДД) = 5-1 /(Д)5, /(Д) = (Ж)). (13)

Таким образом, решение системы (9) запишется следующим образом:

Ф(Н) = 5-

/ехр(гшМ1) 0 0 0

1 0 ехр(гшМ2) 0 0

0 0 ехр(г Лз) 0

0 0 0 ехр(гшЫ4))

5Ф(0). (14)

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

Фт|А,а = Р(Н)(Ф/ + ФД)|А,а , (15)

где Фт — вектор, состоящий их компонент поля, характеризующих прошедшую волну, Фд — вектор, описывающий отражённую волну, а Ф/ — вектор, описывающий падающую волну. Отметим, что каждый из этих векторов состоит из четырёх компонент — двух тангенциальных компонент электрического поля и двух тангенциальных компонент магнитного поля. Р(К) — матрица Берремана для данного образца. Заметим, что для каждой длины волны из необходимого для исследований диапазона и для каждого угла падения задачу составления системы дифференциальных уравнений и нахождения матрицы Берремана надо решать отдельно. Фактически, в (15) при помощи матрицы Берремана были связаны между собой тангенциальные компоненты электромагнитного поля в верхней полуплоскости (над слоем) и в нижней полуплоскости (под слоем). Выражение (15) является системой линейных уравнений. Однако нетрудно заметить, что уравнений там четыре, а неизвестных восемь (по 4 компоненты векторов Фд и Ф^). Чтобы сократить количество неизвестных, необходимо выписать соотношения между тангенциальными компонентами электрического и магнитного полей в верхней и нижней средах [6]:

Hix = -niEiy cosai — -ГуЕу, HIy

HRx = niEiy' cos ai — ryRy,

H

Ry

= niEx cos ai ni Ei

- fX EX ,

— fx Rx

cos ai

(16)

HTx = -n2E2y cos a2 - -r**Ty, HTy = - П2Е2х — r*xTx, y y y y cos a2 ж

rx = ni/ cos ai, ry = ni cos ai, rx = n2/ cos a2, r* = n2 cos a2. Теперь уравнение (15) можно записать в следующем виде:

Фт

( тх \

Ту

' УТУ

V r*xTx J

= Р(h) [Ф/ + Фд] , Ф/ =

\

/ Ex

Ey

—ryEy

у Т'xEx J

Ф

R =

/ Rx

Ry Г yRy

\ TxRx j

(17)

Эта система путём алгебраических преобразований может быть сведена к виду Ах = В, где вектор неизвестных включает в себя компоненты электрических

полей Rr

Ry 1 ТХ 1

Ту. Эту систему можно решать любым устойчивым методом.

Автором используется широко известный численный метод Ьи разложения матрицы коэффициентов системы [3] с дальнейшим вычислением вектора неизвестных по методу Гаусса.

Энергетические коэффициенты отражения (И,) и пропускания (Т) находятся как отношения потоков энергии отражённой электромагнитной волны и прошедшей волны к потоку энергии падающей волны соответственно [7]:

R

4-к

E R, H R с 4тт E т, H т

L -, T = J

с

4тг

Е i, Н i

с

4ж

Е I, Н I

(18)

Вектора с индексом I относятся к падающему свету, с индексом К — к отражённому, а с индексом Т — к прошедшему. Учитывая (16) и правила вычисления векторного произведения, формулы для вычисления энергетических коэффициентов принимают следующий вид:

R

+ IRI2

Ех

cos ai

+ \Ey |2

T

n2 cos a2

+ \Ty I2

ni cosai

Ex

cos ai

+ \Ey |2

(19)

2

2

R

X

X

cos ai

cos a2

2

2

В случае, если коэффициенты преломления для верхней и нижней среды одинаковы ( п1 = п2,а1 = а2), формула для вычисления коэффициента Т принимает следующий вид:

T

+ I Ту |2

Ех

cos ai

+ |Ду I2

(20)

X

cos ai

2

Литература

1. Агранович В. М., Гинзбург В. Л. Кристаллооптика с учетом пространственной дисперсии и теория экситонов. — М.: Наука, 1965. — 376 с.

2. Berreman D. W. Optics in Stratified and Anisotropic Media: 4x4-Matrix Formulation // Journal of the Optical Society of America. — 1972. — No 4. — Pp. 502-510.

3. Уилкинсон Дж., Райнш К. Справочник алгоритмов на языке Алгол. — М.: Машиностроение, 1976. — 389 с.

4. Беллман Р. Введение в теорию матриц. — М.: Мир, 1969. — 368 с.

5. Гантмахер Ф. Р. Теория матриц. — М.: Наука, 1966. — 576 с.

6. Палто С. П. Алгоритм решения оптической задачи для слоистых анизотропных сред // ЖЭТФ. — 2001. — № 4. — С. 638-648.

7. Севастьянов А. Л, Ловецкий К. П., Хохлов А. А. Моделирование прохождения света через многослойные анизотропные покрытия. Восстановление тензора диэлектрической проницаемости для однослойного анизотропного образца // материалы международной научной конференции «Моделирование нелинейных процессов и систем». — 2008. — № 12. — С. 484-498.

UDC 621.39

Modeling of Interaction of Visible Spectra Electromagnetic Wave with Anisotropic Layer

R. P. Lovetskiy, A. A. Khokhlov

Telecommunication Systems Department

Peoples' Friendship University of Russia 6, Miklukho-Maklaya str., 117198, Moscow, Russia

Modeling of interaction of TE and TM light polarizations with optical structures, which characteristic dimensions comparable to the optical wavelength, is actual nowadays because of modern directions in development of microelectronics, optical electronics and medicine. Energy saving solar technologies are developed in USA and Europe very intensively, development and optimization of solar batteries need robust mathematical models of interaction of light and solar cells coverage. These models are also applied in LED and OLED displays development and improvement. High quality mirrors, polarizers and other optical elements (less than 1 micrometer in width) are claimed in modern optical-electronic devices. It is significant that different prosthetic devices can be developed using these models, for example - artificial eye based on multilayer collagen structures.

Key words and phrases: anisotropic optics, Berreman matrix, propagation of electromagnetic wave, boundary conditions for electromagnetic field, Maxwell equations, frequency dispersion.