Моделирование взаимодействия произвольного светового поля с дифракционной решеткой методом Монте-Карло
УДК 535.131: 535.137 + 535.42
МОДЕЛИРОВАНИЕ ВЗАИМОДЕЙСТВИЯ ПРОИЗВОЛЬНОГО СВЕТОВОГО ПОЛЯ С ДИФРАКЦИОННОЙ РЕШЕТКОЙ МЕТОДОМ МОНТЕ-КАРЛО
© 2012 г. В. В. Савуков; И. В. Голубенко, канд. физ.-мат. наук
Балтийский государственный технический университет “Военмех” им. Д.Ф. Устинова, Санкт-Петербург
E-mail: vladimir@savukov.ru
При рассмотрении некоторых специальных вопросов статистической физики возникла необходимость в высокоточном решении задачи дифракции. В настоящей статье сообщается о создании компьютерной программы, служащей инструментальным средством для расчета параметров дифракционных явлений при теоретическом исследовании сложных оптических систем. Программа выполняет решение задачи дифракции строгим методом на основе уравнений Максвелла при заданных граничных условиях. Допускается произвольная, например, диффузная конфигурация исходного светового поля. В качестве дифракционных оптических элементов рассматриваются отражательные решетки с линейным или скрещенным синусоидальным профилем микрорельефа поверхности. При наличии в системе нескольких дифракционных элементов возможен расчет характеристик самосогласованного итогового светового поля.
Ключевые слова: дифракция, поляризация, индикатриса, рассеяние, диффузный.
Коды OCIS: 050.1940, 260.1960, 260.5430, 290.2648, 290.5855
Поступила в редакцию 01.02.2012
Введение
В ряде практически значимых случаев возникает необходимость решения задачи дифракции, представленной в весьма общей форме. Например, исходное световое поле, взаимодействующее с неким дифракционным оптическим элементом (ДОЭ), может иметь любые физически возможные параметры. В частности, такое первичное (до дифракции) излучение может характеризоваться конкретным спектральным составом, а также заданным угловым распределением энергетической яркости и степени поляризации в геометрическом пространстве. Или помещенный в исходное световое поле ДОЭ представляет собой выполненную из назначенного материала отражательную дифракционную решетку с заданной макроскопической формой и параметрами микрорельефа поверхности. В частном случае эта поверхность может быть зеркально гладкой. Если в рассматриваемом пространстве находится более одного ДОЭ, то возникает необходимость в расчете самосогласованного стационарного состояния результирующего светового поля.
Подобные нетривиальные задачи очень редко имеют точный аналитический вариант своего математического описания. Чаще всего бывает, что даже если некий аналитический вариант существует, то он, как правило, получен с применением достаточно произвольных приближений, последующее влияние которых на адекватность создаваемой модели невозможно объективно оценить.
Для решения указанных задач в наибольшей степени подходят методы имитационного моделирования, в частности, метод статистических испытаний, известный как метод Монте-Карло. При этом результат решения весьма сложной исходной задачи сводится к статистическому обобщению результатов решения множества неких более простых (элементарных) задач, каждая из которых разрешима уже непосредственно [1].
Постановка задачи
При рассмотрении некоторых специальных вопросов, связанных с исследованием аксиоматических принципов статистической физики [2–4], возникла необходимость в высокоточном
10 “Оптический журнал”, 79, 7, 2012
исполнении строгого метода решения задачи дифракции, основанного на решении уравнений Максвелла при заданных граничных условиях, который отвечал бы требованиям
1. Метод должен быть применим к задаче взаимодействия произвольного светового поля (в частности, излучения черного тела) с заданной дифракционной решеткой.
2. Элементарный физический объект, дифракция которого подлежит дальнейшему анализу, представляет собой монохромную электромагнитную волну с плоским фронтом распространения. Данный волновой объект может служить квазиклассической моделью так называемого “свободного” фотона, находящегося в зоне, удаленной от источника излучения. В каждом отдельном случае такой единичный фотон характеризуется следующими параметрами:
– длина волны и, как следствие, величина планковской энергии,
– угол падения фотона на поверхность дифракционной решетки (т. е. угол между волновым вектором фотона и нормалью к макроскопической поверхности),
– азимутальный (конический) угол поворота плоскости падения фотона,
– фаза поляризации (сдвиг фазы между электрическими векторами разных поляризаций),
– соотношение между сторонами прямоугольника, описанного вокруг эллипса поляризации.
3. Дифракционная решетка должна представлять собой отражательную решетку фазового типа с одномерным (линейным) или двумерным (скрещенным) синусоидальным микрорельефом поверхности и иметь численные параметры
– период (шаг) решетки – отдельно для каждого измерения,
– полная глубина микрорельефа решетки – отдельно для каждого измерения,
– комплексная диэлектрическая проницаемость реального (т. е. с конечной проводимостью) материала решетки, функционально зависящая от длины волны излучения. Случай идеальной проводимости должен быть описан отдельно.
4. Корректный расчет рассеяния каждого единичного волнового объекта (фотона) на дифракционной решетке должен быть возможен для максимально общего случая сочетаний параметров фотона и решетки. Это, в первую
очередь, относится к анализу явлений в зонах аномалий Рэлея–Вуда, особенно для субволновых решеток.
5. Результаты расчета дифракции единичного фотона должны включать в себя вероятность поглощения данного фотона решеткой и число порядков (когерентных каналов) рассеяния фотона. Для каждого выявленного порядка вычисляются значение эффективности, угловые характеристики отражения от решетки, а также параметры поляризации излучения в данном конкретном канале рассеяния.
Общая методология расчетов
Требуемая структура исходного светового поля (см. п. 1) генерируется с помощью вероятностно обоснованных выборок начальных параметров для каждого отдельного волнового объекта (см. п. 2). Результаты решений (см. пп. 4–5), получаемые для всех элементарных актов дифракционного рассеяния фотонов на решетке (см. п. 3), запоминаются программой. Итоговые параметры излучения находятся в численной форме как статистические моменты распределений, получаемых вследствие обобщения информации расчетов для всех единичных элементарных актов рассеяния. Для большей наглядности итоговая информация выводится также и в форме графических изображений, например, в виде угловых диаграмм яркости рассеянного излучения.
Расчет единичного (элементарного) акта дифракции
Решение поставленной задачи потребовало рассмотрения дифракции на поверхности, которая обладает периодичностью в направлении двух осей декартовой системы координат. Такая поверхность, получившая название трехмерно модулированного профиля или скрещенных дифракционных решеток, образуется при наложении двух двумерно модулированных профилей или обычных дифракционных решеток, каждая из которых периодична вдоль одного из двух взаимно перпендикулярных направлений.
Предполагается проведение численного анализа свойств таких решеток в широком диапазоне изменения отношения длины волны к периоду, включающему резонансную об-
“Оптический журнал”, 79, 7, 2012
11
ласть, в которой это отношение приближается к единице. В резонансной области на характер дифракции в существенной степени оказывает влияние поляризация падающего света. Поэтому для проведения численного анализа может быть пригоден строгий метод расчета, основанный на решении волновых уравнений или же на решении уравнений Максвелла, позволяющий учесть характер поляризации.
Традиционно, строгие теории, на основе которых можно построить вычислительные алгоритмы для расчета дифракционных процессов, делятся на две основные категории: интегральные и дифференциальные. Интегральные методы предполагают численное решение интегральных уравнений или систем связанных интегральных уравнений, а дифференциальные методы приводят к решению бесконечных систем связанных дифференциальных уравнений. Как показала практика применения указанных теорий, при решении систем дифференциальных уравнений возникают ситуации нестабильности, в которых только применение интегральных методов позволяет получить правильные результаты [5]. С другой стороны, возможности применения существующих в настоящее время интегральных методов сильно лимитированы ввиду сложности выполнения граничных условий при произвольном направлении падающих волн. Поэтому выбор метода для проведения данного исследования оказался не простой проблемой.
В настоящее время решению задачи дифракции на поверхности, обладающей двойной периодичностью, уделяется большое внимание и посвящено много публикаций. Есть много различных подходов к решению этой задачи. Существующие строгие методы включают метод конечных разностей, метод вариации границ, анализ на основе связанных волн, а также различные методы использования криволинейных координат, из которых наиболее надежным и эффективным признан С-метод Шандезона [6–8]. Главной особенностью С-метода является использование такого преобразования системы координат, которое превращает поверхность решетки в одну из координатных поверхностей. Данное преобразование называют трансляционной системой координат.
С помощью тензорной формы записи уравнений Максвелла в трансляционной системе координат, которая подробно рассматривается
в работе [9], можно получить линейную систему дифференциальных уравнений с постоянными коэффициентами. Последнее обстоятельство оказывается очень важным, поскольку позволяет находить решение задачи дифракции путем вычисления собственных значений и собственных векторов матрицы этой системы. Данный метод имеет в значительной мере более широкую область применения по сравнению с другими дифференциальными методами и лишен присущих им недостатков. Преобразование рельефной поверхности решетки к формально плоскому виду значительно упрощает запись граничных условий. Метод также позволяет проводить расчеты одновременно для случаев TM и TE поляризации.
Для проведения расчетов на основе С-метода был разработан реализованный на Фортране численный алгоритм, позволяющий проводить достоверные вычисления для различных диапазонов длин волн, параметров профиля и разных видов материала решетки, включая вариант и с неограниченной проводимостью. Прототипом для данного алгоритма послужили работы [10, 11] по развитию С-метода применительно к скрещенным дифракционным решеткам и к решеткам со слоями диэлектрических покрытий.
Кроме дифракционной эффективности, которая является одной из важнейших характеристик решетки, с помощью разработанного алгоритма можно вычислить комплексные амплитуды и фазы для TM и TE поляризованных дифракционных волн, что, в конечном счете, позволяет определить параметры эллиптической поляризации каждой волны, т. е. определить форму и наклон осей эллипса.
Помимо решения задачи для трехмерно модулированной металлической поверхности алгоритм позволяет выполнять расчет рассеянного светового поля в режиме конической дифракции для обычной (линейной) решетки с конечной проводимостью. Кроме того, при незначительной модификации алгоритм может быть использован для прозрачных диэлектрических решеток, а также решеток с диэлектрическим покрытием.
Структура программного комплекса
Созданный программный комплекс состоит из следующих основных компонентов:
– Базовая часть программы, написанная на языке интерпретатора математической систе-
12 “Оптический журнал”, 79, 7, 2012
мы Maple™ версии 15.01 (производство фирмы Waterloo Maple Inc.). Эта часть, имеющая форму интерактивного “рабочего листа” (Worksheet), содержит полную информацию о поставленной задаче: параметры исходного светового поля, все характеристики дифракционных оптических элементов, помещенных в это поле и др. Помимо этого данный рабочий лист имеет исчерпывающее описание того, какую расчетную информацию требуется накапливать, как эта информация должна обрабатываться и в какой визуальной форме следует представлять результаты обработки пользователю.
– Программные модули в форме dllбиблиотек, вызываемые по мере необходимости из базовой части комплекса. Эти модули, исходные тексты которых написаны на языке Fortran 95, содержат алгоритмы расчетов единичных актов дифракции. Они могут использоваться для разных вариантов размерности микрорельефа и проводимости отражающей поверхности ДОЭ. Модули получены с помощью компилятора Compaq Visual Fortran™ версии 6.6 и включают обращения к математической библиотеке IMSL® Fortran Numerical Library версии 4.01.
– Текстовые файлы с таблицами специальных многомерных последовательностей Соболя [12, 13], предназначенные для использования метода Монте-Карло с особо высоким разрешением. Обращение к этим файлам проводится из базовой части комплекса при формировании сетки значений входных параметров задачи. Файлы таблиц имеют неизменный вид для различных вариантов решаемых задач. Поэтому данные файлы были однократно сгенерированы с помощью соответствующей программы на языке Pascal, выполненной в среде системы Borland Delphi™ версии 6.0 фирмы Borland Software Corporation™.
Таким образом, ввод всех исходных данных и вывод получаемых результатов осуществляется через интерактивный рабочий лист (система Maple™) базовой части созданного программного комплекса, функционирующего на платформе Win32®.
Верификация программы
В качестве основного критерия, позволяющего контролировать правильную работу программного алгоритма, было выбрано условие соблюдения энергетического баланса. Для этого сравнивалась энергия, приходящая к
решетке вместе с падающей волной, с энергией, которая уносилась распространяющимися дифракционными волнами или поглощалась в материале решетки. При корректной работе программы энергетический баланс должен всегда выполняться с высокой точностью.
Известно, однако, что выполнение энергетического баланса не может служить единственным доказательством правильности решения задачи дифракции. Поэтому была проведена доскональная проверка разработанного алгоритма, которая заключалась в сравнении результатов расчетов с аналогичной информацией, получаемой интегральными методами как для идеально проводящих решеток, так и для решеток с конечной проводимостью [5]. При этом в числе прочих рассматривались и такие особые режимы, как, например, случай полного поглощения излучения поверхностью решетки [14]. Было также выполнено сравнительное тестирование данных, получаемых для скрещенных решеток, с уже известными результатами, представленными в работе [10].
Для большей уверенности в корректности созданной программы была осуществлена прямая экспериментальная проверка результатов ее работы в одном из самых сложных режимов: тестировался расчет дифракционного рассеяния монохромного диффузного излучения (длина волны = 470 нм) на линейной субволновой фазовой решетке из алюминия (шаг d = 372 нм, полная глубина синусоидального профиля микрорельефа h = 118 нм) при наличии зон интенсивного проявления аномалий Рэлея–Вуда [15, 16].
Согласно описанной логике имитационного моделирования, данное диффузное излучение имеет дискретную структуру (фотонный газ) и определяется как неполяризованное некогерентное электромагнитное излучение, для отдельных фотонов которого с равной вероятностью выполняется любая возможная угловая ориентация в трехмерном геометрическом пространстве их волновых k-векторов. Каждый фотон находится в удаленной от источника зоне и поэтому характеризуется плоским волновым фронтом.
На рис. 1а приведено графическое изображение компьютерной модели пространственной индикатрисы суммарной энергетической яркости рассеянного светового поля. Индикатриса построена в полярной системе координат таким образом, что ее центр соответствует
“Оптический журнал”, 79, 7, 2012
13
(а) (б)
1,5
1
0,5
–1,5 –1 –0,5 0 –0,5 –1 –1,5
0,5 1 1,5
Рис. 1. Сравнительные изображения прогнозируемого и реального световых полей.
нулевому значению угла отражения при обзоре поверхности дифракционной решетки. Величина угла отражения пропорциональна полярному радиусу, и на периферии круговой диаграммы значение этого угла приближается к 90. Азимутальный угол наблюдения поверхности решетки определяется полярным углом диаграммы.
На рис. 1б приведена фотография индикатрисы рассеянного светового поля в реальной физической установке, прогноз работы которой был представлен на рис. 1а. Здесь применен принцип получения трехмерного изображения, фиксирующего параметры объекта сразу почти для всех возможных сочетаний углов наблюдения, известный как “эффект рыбьего глаза”. Данный принцип основан на фотографировании отражения поверхности дифракционной решетки в зеркальном шаре, расположенном на малом расстоянии от этой решетки. Съемка выполняется через небольшое сквозное отверстие в теле самой решетки. Поэтому в центре, справа и слева на фотографии 1б видны три вторичных отражения поверхности решетки от шара, которые не имеют отношения к собственно индикатрисе рассеяния.
Очевидное совпадение общего вида индикатрис на рис. 1а, б говорит об успешном характере проведенной экспериментальной проверки.
Примеры результатов работы
Как уже отмечено ранее, результаты работы описываемого программного комплекса, характеризующие итоговые параметры рас-
сеянного излучения, могут быть представлены в численном виде и в форме графических изображений. Поскольку графический вариант обладает большей наглядностью, то здесь для примера приведены некоторые угловые диаграммы энергетической яркости излучения, отражаемого от металлических поверхностей с разными свойствами.
На рис. 2, 3 представлены графические изображения индикатрис различных компонент энергетической яркости рассеянного светового поля для некоторых видов отражающих поверхностей. Во всех случаях исходное световое поле представляет собой монохроматическое диффузное излучение с = 555 нм и общим числом фотонов в каждом статистическом испытании 218 = 262144. Характер построения этих индикатрис аналогичен тому, который был ранее описан для графика на рис. 1а.
Индикатрисы 2а–в описывают угловое распределение яркости диффузного светового поля, отражаемого от идеально проводящего металлического зеркала (S + P, Sи P-поляризованные компоненты).
Каждый отдельный график автоматически масштабируется так, чтобы максимальным образом выявлять все имеющиеся контрасты плотности рассеянного светового потока. Очевидно, что на изображениях индикатрис 2а–в присутствуют лишь бессистемные проявления флуктуаций рассеянного поля, не образующие каких-либо визуально наблюдаемых макроскопических градиентов. Иначе говоря, при отражении диффузного излучения от идеально проводящего зеркала для всех рассеиваемых компонент выполняется закон Ламберта, что вполне ожидаемо.
Индикатрисы 2г–е описывают угловое распределение яркости диффузного светового поля, отражаемого от медного идеально ровного металлического зеркала. В отличие от графиков 2а–в, на изображениях индикатрис 2г–е присутствуют макроскопические градиенты, масштаб которых для большей наглядности автоматически увеличивается до максимально возможной величины1. При этом флуктуационные явления становятся менее заметны на фоне наблюдаемых макроэффектов.
Индикатрисы 3а–в описывают угловое распределение яркости диффузного светового
1 Их реальный масштаб вычисляется через соответствующие статистические моменты.
14 “Оптический журнал”, 79, 7, 2012
(а)
1,5
1
0,5
(б)
1,5
1
0,5
(в)
1,5
1
0,5
–1,5 –1 –0,5 0 –0,5 –1
0,5
1
11,,55 –1,5
–1
–0,5 0
–0,5
–1
0,5 1
1,5 –1,5
–1
–0,5
0
–0,5
–1
0,5 1 1,5
–1,5
(г)
1,5
–1,5
(д)
1,5
–1,5
(е)
1,5
1 11
0,5 0,5 0,5
–1,5 –1 –0,5 0 –0,5 –1 –1,5
0,5
1
1,5 –1,5
–1
–0,5 0
–0,5
–1
–1,5
0,5
1
1,5 –1,5
–1
–0,5 0
–0,5
–1
–1,5
0,5 1 1,5
Рис. 2. Яркость поверхности идеального (а–в) и медного (г–е) зеркал в диффузном свете. а, г – графики суммарной яркости S- и P-поляризованных отраженных компонент, б, д – графики яркости S-поляризованной составляющей отраженного излучения, в, е – графики яркости P-поляризованной составляющей отраженного излучения.
(а)
1,5
1
0,5
(б)
1,5
1
0,5
(в)
1,5
1
0,5
–1,5 –1
–0,5 0
–0,5
0,5 1 1,5 –1,5 –1 –0,5 0 –0,5
0,5
1
1,5 –1,5
–1 –0,5 0
–0,5
0,5 1 1,5
–1
–1,5
(г)
1,5
1
–1
–1,5
(д)
1,5
1
–1
–1,5
(е)
1,5
1
0,5 0,5
0,5
–1,5 –1 –0,5 0 –0,5
0,5
1
1,5 –1,5
–1 –0,5
0
–0,5
0,5 1 1,5 –1,5 –1 –0,5 0 –0,5
0,5 1 1,5
–1 –1
–1
–1,5 –1,5 –1,5
Рис. 3. Сравнение эффективности методов генерирования сетки исходных данных. Метод последовательностей Соболя (а–в), стандартный генератор псевдослучайных чисел, входящий в состав математической системы Maple™ (г–е). а, г – графики суммарной яркости S- и P-поляризованных отраженных компонент, б, д – графики яркости S-поляризованной составляющей отраженного излучения, в, е – графики яркости P-поляризованной составляющей отраженного излучения.
“Оптический журнал”, 79, 7, 2012
15
поля, дифракционно рассеиваемого на линейной фазовой решетке из меди (d = 544 нм, h = = 337 нм).
Индикатрисы 3г–е описывают те же самые угловые распределения яркости, что и представленные на рис. 3а–в. Отличие состоит лишь в том, что в случае 3а–в для формирования значений сетки исходных данных статистических испытаний использовался высокоэффективный метод последовательностей Соболя, а варианты 3г–е были получены с использованием для указанных целей стандартного генератора псевдослучайных чисел, входящего в состав математической системы Maple™.
Здесь не приводятся результаты расчетов, полученные для сетки с равномерным шагом разбиения области допустимых значений входных параметров, поскольку данный тип разбиения является, как известно, наихудшим для многомерных областей.
Выводы
Широкие функциональные возможности созданного программного комплекса, а также принципиально высокая достоверность получаемых с его помощью результатов, позволяют эффективно использовать данный инструмент в теоретических исследованиях, характерных для физической оптики. Более того, во многих случаях, связанных и с экспериментальным изучением нетривиальных оптических систем, становится возможной замена прямого физического эксперимента на имитационный, реализуемый на компьютере. Чаще всего такая замена оправдана с экономической точки зрения. Кроме того, она дает выигрыш во времени проведения работ.
Авторы выражают глубокую благодарность доктору Даниэлю Мэстру (Франция) за конструктивные замечания, способствовавшие выполнению данного проекта.
*****
ЛИТЕРАТУРА
1. Соболь И.М. Численные методы Монте-Карло. М.: Наука, 1973. 311 с.
2. Савуков В.В. Нарушение изотропности диффузного излучения вследствие его дифракции на многомерных регулярных структурах // Оптический журнал. 2010. Т. 77. № 1. С. 95–100.
3. Савуков В.В. Уточнение аксиоматических принципов статистической физики // Деп. в ВИНИТИ. № 1249B2004 от 16.07.2004. URL: http://www.savukov.ru/viniti_1249_b2004_full_rus.pdf
4. Савуков В.В. Нарушение закона Ламберта при дифракции диффузного фотонного газа на многомерных регулярных структурах // Деп. в ВИНИТИ. № 507-B2009 от 03.08.2009. URL: http://www.savukov.ru/ viniti_0507_b2009_full_rus.pdf
5. Electromagnetic Theory of Gratings / Ed. by Petit R. Berlin, N. Y.: Springer-Verlag, 1980. 284 p.
6. Chandezon J., Cornet G., Raoult G. Propagation des ondes dans les guides cylindriques à génératrices sinusoïdales // C. R. Acad. Sci. (Paris). 1973. № 276B. P. 507–509.
7. Chandezon J., Maystre D., Raoult G. A new theoretical method for diffraction gratings and its numerical application // J. Optics (Paris). 1980. V. 11. № 4. P. 235–241.
8. Chandezon J., Dupuis M.T., Cornet G., Maystre D. Multicoated gratings: a differential formalism applicable in the entire optical region // J. Opt. Soc. Am. 1982. V. 72. № 7. P. 839–847.
9. Post E. J. Formal structure of electromagnetics: General covariance and electromagnetics. Amsterdam: NorthHolland Pub. Co., N. Y.: Interscience Publishers, 1962. Series in physics. 224 p.
10. Granet G. Diffraction par des surfaces biperiodiques: resolution en coordonnees non orthogonales // Pure Appl. Opt. 1995. V. 4. P. 777–793.
11. Granet G. Analysis of diffraction by surface-relief crossed gratings with use of the Chandezon method: application to multilayer crossed gratings // J. Opt. Soc. Am. 1998. V. 15. № 5. P. 1121–1131.
12. Соболь И.М. Многомерные квадратурные формулы и функции Хаара. М.: Наука, 1969. 288 с.
16 “Оптический журнал”, 79, 7, 2012
13. Соболь И.М., Статников Р.Б. Выбор оптимальных параметров в задачах со многими критериями. М.: Наука, 1981. 110 с.
14. Maystre D., Petit R. Brewster incidence for metallic gratings // Opt. Commun. 1976. V. 17. № 2. P. 196–200.
15. Wood R.W. On a remarkable case of uneven distribution of light in a diffraction grating spectrum // Philosophical Magazine. 1902. V. 4. P. 396–402.
16. Lord Rayleigh. On the Dynamical Theory of Gratings // Proc. Royal Soc. London. 1907. Series A 79. P. 399–416.
“Оптический журнал”, 79, 7, 2012
17
МОДЕЛИРОВАНИЕ ВЗАИМОДЕЙСТВИЯ ПРОИЗВОЛЬНОГО СВЕТОВОГО ПОЛЯ С ДИФРАКЦИОННОЙ РЕШЕТКОЙ МЕТОДОМ МОНТЕ-КАРЛО
© 2012 г. В. В. Савуков; И. В. Голубенко, канд. физ.-мат. наук
Балтийский государственный технический университет “Военмех” им. Д.Ф. Устинова, Санкт-Петербург
E-mail: vladimir@savukov.ru
При рассмотрении некоторых специальных вопросов статистической физики возникла необходимость в высокоточном решении задачи дифракции. В настоящей статье сообщается о создании компьютерной программы, служащей инструментальным средством для расчета параметров дифракционных явлений при теоретическом исследовании сложных оптических систем. Программа выполняет решение задачи дифракции строгим методом на основе уравнений Максвелла при заданных граничных условиях. Допускается произвольная, например, диффузная конфигурация исходного светового поля. В качестве дифракционных оптических элементов рассматриваются отражательные решетки с линейным или скрещенным синусоидальным профилем микрорельефа поверхности. При наличии в системе нескольких дифракционных элементов возможен расчет характеристик самосогласованного итогового светового поля.
Ключевые слова: дифракция, поляризация, индикатриса, рассеяние, диффузный.
Коды OCIS: 050.1940, 260.1960, 260.5430, 290.2648, 290.5855
Поступила в редакцию 01.02.2012
Введение
В ряде практически значимых случаев возникает необходимость решения задачи дифракции, представленной в весьма общей форме. Например, исходное световое поле, взаимодействующее с неким дифракционным оптическим элементом (ДОЭ), может иметь любые физически возможные параметры. В частности, такое первичное (до дифракции) излучение может характеризоваться конкретным спектральным составом, а также заданным угловым распределением энергетической яркости и степени поляризации в геометрическом пространстве. Или помещенный в исходное световое поле ДОЭ представляет собой выполненную из назначенного материала отражательную дифракционную решетку с заданной макроскопической формой и параметрами микрорельефа поверхности. В частном случае эта поверхность может быть зеркально гладкой. Если в рассматриваемом пространстве находится более одного ДОЭ, то возникает необходимость в расчете самосогласованного стационарного состояния результирующего светового поля.
Подобные нетривиальные задачи очень редко имеют точный аналитический вариант своего математического описания. Чаще всего бывает, что даже если некий аналитический вариант существует, то он, как правило, получен с применением достаточно произвольных приближений, последующее влияние которых на адекватность создаваемой модели невозможно объективно оценить.
Для решения указанных задач в наибольшей степени подходят методы имитационного моделирования, в частности, метод статистических испытаний, известный как метод Монте-Карло. При этом результат решения весьма сложной исходной задачи сводится к статистическому обобщению результатов решения множества неких более простых (элементарных) задач, каждая из которых разрешима уже непосредственно [1].
Постановка задачи
При рассмотрении некоторых специальных вопросов, связанных с исследованием аксиоматических принципов статистической физики [2–4], возникла необходимость в высокоточном
10 “Оптический журнал”, 79, 7, 2012
исполнении строгого метода решения задачи дифракции, основанного на решении уравнений Максвелла при заданных граничных условиях, который отвечал бы требованиям
1. Метод должен быть применим к задаче взаимодействия произвольного светового поля (в частности, излучения черного тела) с заданной дифракционной решеткой.
2. Элементарный физический объект, дифракция которого подлежит дальнейшему анализу, представляет собой монохромную электромагнитную волну с плоским фронтом распространения. Данный волновой объект может служить квазиклассической моделью так называемого “свободного” фотона, находящегося в зоне, удаленной от источника излучения. В каждом отдельном случае такой единичный фотон характеризуется следующими параметрами:
– длина волны и, как следствие, величина планковской энергии,
– угол падения фотона на поверхность дифракционной решетки (т. е. угол между волновым вектором фотона и нормалью к макроскопической поверхности),
– азимутальный (конический) угол поворота плоскости падения фотона,
– фаза поляризации (сдвиг фазы между электрическими векторами разных поляризаций),
– соотношение между сторонами прямоугольника, описанного вокруг эллипса поляризации.
3. Дифракционная решетка должна представлять собой отражательную решетку фазового типа с одномерным (линейным) или двумерным (скрещенным) синусоидальным микрорельефом поверхности и иметь численные параметры
– период (шаг) решетки – отдельно для каждого измерения,
– полная глубина микрорельефа решетки – отдельно для каждого измерения,
– комплексная диэлектрическая проницаемость реального (т. е. с конечной проводимостью) материала решетки, функционально зависящая от длины волны излучения. Случай идеальной проводимости должен быть описан отдельно.
4. Корректный расчет рассеяния каждого единичного волнового объекта (фотона) на дифракционной решетке должен быть возможен для максимально общего случая сочетаний параметров фотона и решетки. Это, в первую
очередь, относится к анализу явлений в зонах аномалий Рэлея–Вуда, особенно для субволновых решеток.
5. Результаты расчета дифракции единичного фотона должны включать в себя вероятность поглощения данного фотона решеткой и число порядков (когерентных каналов) рассеяния фотона. Для каждого выявленного порядка вычисляются значение эффективности, угловые характеристики отражения от решетки, а также параметры поляризации излучения в данном конкретном канале рассеяния.
Общая методология расчетов
Требуемая структура исходного светового поля (см. п. 1) генерируется с помощью вероятностно обоснованных выборок начальных параметров для каждого отдельного волнового объекта (см. п. 2). Результаты решений (см. пп. 4–5), получаемые для всех элементарных актов дифракционного рассеяния фотонов на решетке (см. п. 3), запоминаются программой. Итоговые параметры излучения находятся в численной форме как статистические моменты распределений, получаемых вследствие обобщения информации расчетов для всех единичных элементарных актов рассеяния. Для большей наглядности итоговая информация выводится также и в форме графических изображений, например, в виде угловых диаграмм яркости рассеянного излучения.
Расчет единичного (элементарного) акта дифракции
Решение поставленной задачи потребовало рассмотрения дифракции на поверхности, которая обладает периодичностью в направлении двух осей декартовой системы координат. Такая поверхность, получившая название трехмерно модулированного профиля или скрещенных дифракционных решеток, образуется при наложении двух двумерно модулированных профилей или обычных дифракционных решеток, каждая из которых периодична вдоль одного из двух взаимно перпендикулярных направлений.
Предполагается проведение численного анализа свойств таких решеток в широком диапазоне изменения отношения длины волны к периоду, включающему резонансную об-
“Оптический журнал”, 79, 7, 2012
11
ласть, в которой это отношение приближается к единице. В резонансной области на характер дифракции в существенной степени оказывает влияние поляризация падающего света. Поэтому для проведения численного анализа может быть пригоден строгий метод расчета, основанный на решении волновых уравнений или же на решении уравнений Максвелла, позволяющий учесть характер поляризации.
Традиционно, строгие теории, на основе которых можно построить вычислительные алгоритмы для расчета дифракционных процессов, делятся на две основные категории: интегральные и дифференциальные. Интегральные методы предполагают численное решение интегральных уравнений или систем связанных интегральных уравнений, а дифференциальные методы приводят к решению бесконечных систем связанных дифференциальных уравнений. Как показала практика применения указанных теорий, при решении систем дифференциальных уравнений возникают ситуации нестабильности, в которых только применение интегральных методов позволяет получить правильные результаты [5]. С другой стороны, возможности применения существующих в настоящее время интегральных методов сильно лимитированы ввиду сложности выполнения граничных условий при произвольном направлении падающих волн. Поэтому выбор метода для проведения данного исследования оказался не простой проблемой.
В настоящее время решению задачи дифракции на поверхности, обладающей двойной периодичностью, уделяется большое внимание и посвящено много публикаций. Есть много различных подходов к решению этой задачи. Существующие строгие методы включают метод конечных разностей, метод вариации границ, анализ на основе связанных волн, а также различные методы использования криволинейных координат, из которых наиболее надежным и эффективным признан С-метод Шандезона [6–8]. Главной особенностью С-метода является использование такого преобразования системы координат, которое превращает поверхность решетки в одну из координатных поверхностей. Данное преобразование называют трансляционной системой координат.
С помощью тензорной формы записи уравнений Максвелла в трансляционной системе координат, которая подробно рассматривается
в работе [9], можно получить линейную систему дифференциальных уравнений с постоянными коэффициентами. Последнее обстоятельство оказывается очень важным, поскольку позволяет находить решение задачи дифракции путем вычисления собственных значений и собственных векторов матрицы этой системы. Данный метод имеет в значительной мере более широкую область применения по сравнению с другими дифференциальными методами и лишен присущих им недостатков. Преобразование рельефной поверхности решетки к формально плоскому виду значительно упрощает запись граничных условий. Метод также позволяет проводить расчеты одновременно для случаев TM и TE поляризации.
Для проведения расчетов на основе С-метода был разработан реализованный на Фортране численный алгоритм, позволяющий проводить достоверные вычисления для различных диапазонов длин волн, параметров профиля и разных видов материала решетки, включая вариант и с неограниченной проводимостью. Прототипом для данного алгоритма послужили работы [10, 11] по развитию С-метода применительно к скрещенным дифракционным решеткам и к решеткам со слоями диэлектрических покрытий.
Кроме дифракционной эффективности, которая является одной из важнейших характеристик решетки, с помощью разработанного алгоритма можно вычислить комплексные амплитуды и фазы для TM и TE поляризованных дифракционных волн, что, в конечном счете, позволяет определить параметры эллиптической поляризации каждой волны, т. е. определить форму и наклон осей эллипса.
Помимо решения задачи для трехмерно модулированной металлической поверхности алгоритм позволяет выполнять расчет рассеянного светового поля в режиме конической дифракции для обычной (линейной) решетки с конечной проводимостью. Кроме того, при незначительной модификации алгоритм может быть использован для прозрачных диэлектрических решеток, а также решеток с диэлектрическим покрытием.
Структура программного комплекса
Созданный программный комплекс состоит из следующих основных компонентов:
– Базовая часть программы, написанная на языке интерпретатора математической систе-
12 “Оптический журнал”, 79, 7, 2012
мы Maple™ версии 15.01 (производство фирмы Waterloo Maple Inc.). Эта часть, имеющая форму интерактивного “рабочего листа” (Worksheet), содержит полную информацию о поставленной задаче: параметры исходного светового поля, все характеристики дифракционных оптических элементов, помещенных в это поле и др. Помимо этого данный рабочий лист имеет исчерпывающее описание того, какую расчетную информацию требуется накапливать, как эта информация должна обрабатываться и в какой визуальной форме следует представлять результаты обработки пользователю.
– Программные модули в форме dllбиблиотек, вызываемые по мере необходимости из базовой части комплекса. Эти модули, исходные тексты которых написаны на языке Fortran 95, содержат алгоритмы расчетов единичных актов дифракции. Они могут использоваться для разных вариантов размерности микрорельефа и проводимости отражающей поверхности ДОЭ. Модули получены с помощью компилятора Compaq Visual Fortran™ версии 6.6 и включают обращения к математической библиотеке IMSL® Fortran Numerical Library версии 4.01.
– Текстовые файлы с таблицами специальных многомерных последовательностей Соболя [12, 13], предназначенные для использования метода Монте-Карло с особо высоким разрешением. Обращение к этим файлам проводится из базовой части комплекса при формировании сетки значений входных параметров задачи. Файлы таблиц имеют неизменный вид для различных вариантов решаемых задач. Поэтому данные файлы были однократно сгенерированы с помощью соответствующей программы на языке Pascal, выполненной в среде системы Borland Delphi™ версии 6.0 фирмы Borland Software Corporation™.
Таким образом, ввод всех исходных данных и вывод получаемых результатов осуществляется через интерактивный рабочий лист (система Maple™) базовой части созданного программного комплекса, функционирующего на платформе Win32®.
Верификация программы
В качестве основного критерия, позволяющего контролировать правильную работу программного алгоритма, было выбрано условие соблюдения энергетического баланса. Для этого сравнивалась энергия, приходящая к
решетке вместе с падающей волной, с энергией, которая уносилась распространяющимися дифракционными волнами или поглощалась в материале решетки. При корректной работе программы энергетический баланс должен всегда выполняться с высокой точностью.
Известно, однако, что выполнение энергетического баланса не может служить единственным доказательством правильности решения задачи дифракции. Поэтому была проведена доскональная проверка разработанного алгоритма, которая заключалась в сравнении результатов расчетов с аналогичной информацией, получаемой интегральными методами как для идеально проводящих решеток, так и для решеток с конечной проводимостью [5]. При этом в числе прочих рассматривались и такие особые режимы, как, например, случай полного поглощения излучения поверхностью решетки [14]. Было также выполнено сравнительное тестирование данных, получаемых для скрещенных решеток, с уже известными результатами, представленными в работе [10].
Для большей уверенности в корректности созданной программы была осуществлена прямая экспериментальная проверка результатов ее работы в одном из самых сложных режимов: тестировался расчет дифракционного рассеяния монохромного диффузного излучения (длина волны = 470 нм) на линейной субволновой фазовой решетке из алюминия (шаг d = 372 нм, полная глубина синусоидального профиля микрорельефа h = 118 нм) при наличии зон интенсивного проявления аномалий Рэлея–Вуда [15, 16].
Согласно описанной логике имитационного моделирования, данное диффузное излучение имеет дискретную структуру (фотонный газ) и определяется как неполяризованное некогерентное электромагнитное излучение, для отдельных фотонов которого с равной вероятностью выполняется любая возможная угловая ориентация в трехмерном геометрическом пространстве их волновых k-векторов. Каждый фотон находится в удаленной от источника зоне и поэтому характеризуется плоским волновым фронтом.
На рис. 1а приведено графическое изображение компьютерной модели пространственной индикатрисы суммарной энергетической яркости рассеянного светового поля. Индикатриса построена в полярной системе координат таким образом, что ее центр соответствует
“Оптический журнал”, 79, 7, 2012
13
(а) (б)
1,5
1
0,5
–1,5 –1 –0,5 0 –0,5 –1 –1,5
0,5 1 1,5
Рис. 1. Сравнительные изображения прогнозируемого и реального световых полей.
нулевому значению угла отражения при обзоре поверхности дифракционной решетки. Величина угла отражения пропорциональна полярному радиусу, и на периферии круговой диаграммы значение этого угла приближается к 90. Азимутальный угол наблюдения поверхности решетки определяется полярным углом диаграммы.
На рис. 1б приведена фотография индикатрисы рассеянного светового поля в реальной физической установке, прогноз работы которой был представлен на рис. 1а. Здесь применен принцип получения трехмерного изображения, фиксирующего параметры объекта сразу почти для всех возможных сочетаний углов наблюдения, известный как “эффект рыбьего глаза”. Данный принцип основан на фотографировании отражения поверхности дифракционной решетки в зеркальном шаре, расположенном на малом расстоянии от этой решетки. Съемка выполняется через небольшое сквозное отверстие в теле самой решетки. Поэтому в центре, справа и слева на фотографии 1б видны три вторичных отражения поверхности решетки от шара, которые не имеют отношения к собственно индикатрисе рассеяния.
Очевидное совпадение общего вида индикатрис на рис. 1а, б говорит об успешном характере проведенной экспериментальной проверки.
Примеры результатов работы
Как уже отмечено ранее, результаты работы описываемого программного комплекса, характеризующие итоговые параметры рас-
сеянного излучения, могут быть представлены в численном виде и в форме графических изображений. Поскольку графический вариант обладает большей наглядностью, то здесь для примера приведены некоторые угловые диаграммы энергетической яркости излучения, отражаемого от металлических поверхностей с разными свойствами.
На рис. 2, 3 представлены графические изображения индикатрис различных компонент энергетической яркости рассеянного светового поля для некоторых видов отражающих поверхностей. Во всех случаях исходное световое поле представляет собой монохроматическое диффузное излучение с = 555 нм и общим числом фотонов в каждом статистическом испытании 218 = 262144. Характер построения этих индикатрис аналогичен тому, который был ранее описан для графика на рис. 1а.
Индикатрисы 2а–в описывают угловое распределение яркости диффузного светового поля, отражаемого от идеально проводящего металлического зеркала (S + P, Sи P-поляризованные компоненты).
Каждый отдельный график автоматически масштабируется так, чтобы максимальным образом выявлять все имеющиеся контрасты плотности рассеянного светового потока. Очевидно, что на изображениях индикатрис 2а–в присутствуют лишь бессистемные проявления флуктуаций рассеянного поля, не образующие каких-либо визуально наблюдаемых макроскопических градиентов. Иначе говоря, при отражении диффузного излучения от идеально проводящего зеркала для всех рассеиваемых компонент выполняется закон Ламберта, что вполне ожидаемо.
Индикатрисы 2г–е описывают угловое распределение яркости диффузного светового поля, отражаемого от медного идеально ровного металлического зеркала. В отличие от графиков 2а–в, на изображениях индикатрис 2г–е присутствуют макроскопические градиенты, масштаб которых для большей наглядности автоматически увеличивается до максимально возможной величины1. При этом флуктуационные явления становятся менее заметны на фоне наблюдаемых макроэффектов.
Индикатрисы 3а–в описывают угловое распределение яркости диффузного светового
1 Их реальный масштаб вычисляется через соответствующие статистические моменты.
14 “Оптический журнал”, 79, 7, 2012
(а)
1,5
1
0,5
(б)
1,5
1
0,5
(в)
1,5
1
0,5
–1,5 –1 –0,5 0 –0,5 –1
0,5
1
11,,55 –1,5
–1
–0,5 0
–0,5
–1
0,5 1
1,5 –1,5
–1
–0,5
0
–0,5
–1
0,5 1 1,5
–1,5
(г)
1,5
–1,5
(д)
1,5
–1,5
(е)
1,5
1 11
0,5 0,5 0,5
–1,5 –1 –0,5 0 –0,5 –1 –1,5
0,5
1
1,5 –1,5
–1
–0,5 0
–0,5
–1
–1,5
0,5
1
1,5 –1,5
–1
–0,5 0
–0,5
–1
–1,5
0,5 1 1,5
Рис. 2. Яркость поверхности идеального (а–в) и медного (г–е) зеркал в диффузном свете. а, г – графики суммарной яркости S- и P-поляризованных отраженных компонент, б, д – графики яркости S-поляризованной составляющей отраженного излучения, в, е – графики яркости P-поляризованной составляющей отраженного излучения.
(а)
1,5
1
0,5
(б)
1,5
1
0,5
(в)
1,5
1
0,5
–1,5 –1
–0,5 0
–0,5
0,5 1 1,5 –1,5 –1 –0,5 0 –0,5
0,5
1
1,5 –1,5
–1 –0,5 0
–0,5
0,5 1 1,5
–1
–1,5
(г)
1,5
1
–1
–1,5
(д)
1,5
1
–1
–1,5
(е)
1,5
1
0,5 0,5
0,5
–1,5 –1 –0,5 0 –0,5
0,5
1
1,5 –1,5
–1 –0,5
0
–0,5
0,5 1 1,5 –1,5 –1 –0,5 0 –0,5
0,5 1 1,5
–1 –1
–1
–1,5 –1,5 –1,5
Рис. 3. Сравнение эффективности методов генерирования сетки исходных данных. Метод последовательностей Соболя (а–в), стандартный генератор псевдослучайных чисел, входящий в состав математической системы Maple™ (г–е). а, г – графики суммарной яркости S- и P-поляризованных отраженных компонент, б, д – графики яркости S-поляризованной составляющей отраженного излучения, в, е – графики яркости P-поляризованной составляющей отраженного излучения.
“Оптический журнал”, 79, 7, 2012
15
поля, дифракционно рассеиваемого на линейной фазовой решетке из меди (d = 544 нм, h = = 337 нм).
Индикатрисы 3г–е описывают те же самые угловые распределения яркости, что и представленные на рис. 3а–в. Отличие состоит лишь в том, что в случае 3а–в для формирования значений сетки исходных данных статистических испытаний использовался высокоэффективный метод последовательностей Соболя, а варианты 3г–е были получены с использованием для указанных целей стандартного генератора псевдослучайных чисел, входящего в состав математической системы Maple™.
Здесь не приводятся результаты расчетов, полученные для сетки с равномерным шагом разбиения области допустимых значений входных параметров, поскольку данный тип разбиения является, как известно, наихудшим для многомерных областей.
Выводы
Широкие функциональные возможности созданного программного комплекса, а также принципиально высокая достоверность получаемых с его помощью результатов, позволяют эффективно использовать данный инструмент в теоретических исследованиях, характерных для физической оптики. Более того, во многих случаях, связанных и с экспериментальным изучением нетривиальных оптических систем, становится возможной замена прямого физического эксперимента на имитационный, реализуемый на компьютере. Чаще всего такая замена оправдана с экономической точки зрения. Кроме того, она дает выигрыш во времени проведения работ.
Авторы выражают глубокую благодарность доктору Даниэлю Мэстру (Франция) за конструктивные замечания, способствовавшие выполнению данного проекта.
*****
ЛИТЕРАТУРА
1. Соболь И.М. Численные методы Монте-Карло. М.: Наука, 1973. 311 с.
2. Савуков В.В. Нарушение изотропности диффузного излучения вследствие его дифракции на многомерных регулярных структурах // Оптический журнал. 2010. Т. 77. № 1. С. 95–100.
3. Савуков В.В. Уточнение аксиоматических принципов статистической физики // Деп. в ВИНИТИ. № 1249B2004 от 16.07.2004. URL: http://www.savukov.ru/viniti_1249_b2004_full_rus.pdf
4. Савуков В.В. Нарушение закона Ламберта при дифракции диффузного фотонного газа на многомерных регулярных структурах // Деп. в ВИНИТИ. № 507-B2009 от 03.08.2009. URL: http://www.savukov.ru/ viniti_0507_b2009_full_rus.pdf
5. Electromagnetic Theory of Gratings / Ed. by Petit R. Berlin, N. Y.: Springer-Verlag, 1980. 284 p.
6. Chandezon J., Cornet G., Raoult G. Propagation des ondes dans les guides cylindriques à génératrices sinusoïdales // C. R. Acad. Sci. (Paris). 1973. № 276B. P. 507–509.
7. Chandezon J., Maystre D., Raoult G. A new theoretical method for diffraction gratings and its numerical application // J. Optics (Paris). 1980. V. 11. № 4. P. 235–241.
8. Chandezon J., Dupuis M.T., Cornet G., Maystre D. Multicoated gratings: a differential formalism applicable in the entire optical region // J. Opt. Soc. Am. 1982. V. 72. № 7. P. 839–847.
9. Post E. J. Formal structure of electromagnetics: General covariance and electromagnetics. Amsterdam: NorthHolland Pub. Co., N. Y.: Interscience Publishers, 1962. Series in physics. 224 p.
10. Granet G. Diffraction par des surfaces biperiodiques: resolution en coordonnees non orthogonales // Pure Appl. Opt. 1995. V. 4. P. 777–793.
11. Granet G. Analysis of diffraction by surface-relief crossed gratings with use of the Chandezon method: application to multilayer crossed gratings // J. Opt. Soc. Am. 1998. V. 15. № 5. P. 1121–1131.
12. Соболь И.М. Многомерные квадратурные формулы и функции Хаара. М.: Наука, 1969. 288 с.
16 “Оптический журнал”, 79, 7, 2012
13. Соболь И.М., Статников Р.Б. Выбор оптимальных параметров в задачах со многими критериями. М.: Наука, 1981. 110 с.
14. Maystre D., Petit R. Brewster incidence for metallic gratings // Opt. Commun. 1976. V. 17. № 2. P. 196–200.
15. Wood R.W. On a remarkable case of uneven distribution of light in a diffraction grating spectrum // Philosophical Magazine. 1902. V. 4. P. 396–402.
16. Lord Rayleigh. On the Dynamical Theory of Gratings // Proc. Royal Soc. London. 1907. Series A 79. P. 399–416.
“Оптический журнал”, 79, 7, 2012
17