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

УДК 541.144

РАСЧЕТ ГИПЕРГЕОМЕТРИЧЕСКИХ МОД

© 2007 С.А. Балалаев1, С.Н. Хонина2, В.В. Котляр2

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

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

Введение

В данной работе рассмотрен один из типов модовых решений параксиального волнового уравнения (или уравнения Шредингера) в цилиндрических координатах. Эти решения описывают семейство световых полей, которые сохраняют свою структуру, изменяясь только масштабно. Так как функции, описывающие эти моды, содержат вырожденную гипергеометрическую функцию, то их называют гипергеометрическими (ГГ) модами. Ранее были исследованы такие решения параксиального волнового уравнения: в декартовой системе координат - моды Эрмита-Гаусса [1], в цилиндрической системе координат - моды Лагерра-Гаусса [2] и моды Бесселя [3], а также в эллиптической системе координат - пучки Айнса-Гаусса, названные так, поскольку их решение получается в виде произведения гауссовой функции на многочлены Айнса [4-8]. Известны также моды Эрмита-Лагерра-Гаусса [9], которые в зависимости от некоторого параметра могут переходить в моды Эрмита-Гаусса или Лагер-ра-Гаусса.

ГГ моды как и параксиальные моды Бесселя обладают бесконечной энергией. На практике их можно сформировать только приближенно и на конечном расстоянии с помощью спиральной фазовой пластины [10-12], освещенной плоской волной. У любой ГГ моды при п£-0 всегда (кроме начальной плоскости)

в центре (на оптической оси) имеется нуль интенсивности, характерный для оптических вихрей [13]. Кроме этого амплитуда ГГ моды в начальной плоскости (г=0) имеет гиперболическую особенность 1/г при г=0.

Подробно гипергеометрические моды рассмотрены в [14, 15]. Однако полученные аналитические решения оказались довольно сложными для программной реализации в связи с численной неустойчивостью прямого расчета гипергеометрических функции, входящих в состав мод. При больших значениях аргумента происходит переполнение значения функции. Основной задачей данной работы является разработка и реализация новых методов для расчета ГГ мод. Также проведено сравнение распространения точных бесконечных ГГ мод в свободном пространстве с модами, ограниченными апертурой и рассчитанных с использованием интеграла Кирхгофа.

Гипергеометрические моды

Параксиальное волновое уравнение в цилиндрических координатах (уравнение типа Шредингера) имеет вид:

д д2 1 5 1

2± — + —- +--+ —

дг дг г дг

2

г2 дф2

Е(г, ф, г) = 0, (1)

где (г j) - поперечные полярные координаты, г- координата, направленная вдоль оптической оси, k=2p//-волновое число света с

длинной волны I. Уравнению (1) удовлетворяют функции, образующие ортонормирован-ный базис:

1

^2) = ~7жп\I кМ2

iy -1 22

ехр

Ш < 14

- — (п -7 +1)

22

п + ¡7 +1

2

, (_)

( п +1_¡7 ¡кг ^

1 —п+1; — I ехрои«^,

где п>0 - целое число, g - комплексное число; данные числа в дальнейшем будем называть параметрами ГГ мод; м - вещественный параметр, задающий масштаб ГГ моды, аналогичен радиусу перетяжки гауссового пучка (в данной статье будут рассмотрены пучки для которых м=1); Г(х) - гамма функция; ^(а, Ь;х) - вырожденная (или конфлюэнтная) гипергеометрическая функция:

1 F1(a, Ъ, х) =

Г(Ъ)

Г(а)Г(Ъ _ а)

| Г-1(1 _ t )Ъ-а-1 ехр( ха И Л(3)

Выражение (3) может также быть записано в виде ряда Тейлора (функция Куммера):

! Fl(a, Ъ, х) = £

(а)тхт

(4)

т=0 (Ъ)тШ\

где (а)т=а(а+ 1)(а+2) .. .(а+т-1)- символ Пох-гаммера, (а)0=1.

В выражениях (2) и (3) используется гамма функция для вещественного аргумента [16]:

Г(х) = I tx_1-~tdt.

(г^2) = I

ехр

¡я., . 1Ч

--(п - ¡7 +1)

4

кг__ 2 2

п + ¡7 +1

2

(

(6)

1 Fl

п\ +1 -17 . . ¡кг

—-, п +1;-

2 11 22

ехр(/'п<р).

Еу„ (г,Ф, 2 = 0) =

ехр(/п<р) ( г

2я

¡7 -1

(7)

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

2 ^кг ( 1 А

Е(х, у, 2) = -2- IIЕ(£, п) -~г (¡к - -у^п , (8)

где г = ^1 (х-|)2 + (у-п)2 + 22 , Е0(х,И) - входное световое поле в декартовых координатах, вычисляемое с помощью (7). Далее выходное поле Е(х,у,2), полученное с использованием (8) сравнивается с полем, вычисленным при помощи аналитических решений (2) и (6).

Алгоритм расчета

Первым шагом в написании алгоритма, реализующего решение (2) было исследование гамма-функции (5). Для данного семейства гамма-функций имеется ряд численных методов для вычисления ее значений на компьютере, однако все они приводят к результату гораздо худшему, чем аппроксимация методом Ланцоша (метод, полученный математиком, занимавшимся симметричной проблемой собственных значений) [17]:

Г(х +1) = (х + и + })х+2 ехр[-(х + и + })]>

х +1 х + 2

- +... + -

х + Ж

- + -

(9)

(5)

Решение (2) можно записать для случая с отрицательными значениями параметра моды п<0:

¡'7-1

(-1)1( 22 ГТ

где N и и - произвольные параметры уравнения, с .ск - весовые коэффициенты, - - погрешность вычислений. Выражение (9) справедливо для х > 0, для других значений аргумента используют формулу отражения [18]:

Г(1 - х) = -

ях

При этом на нулевом расстоянии вдоль оптической оси, при 2 = 0:

Г(х^ш(жх) Г(1 + х^ш(жх) ' (10)

Основной проблемой при исследовании численного решения (9) был подбор оптимальных параметров N и и, при которых можно было вычислить коэффициенты с ,ск. Для этого была построена система линейных алгебраических уравнений с N неизвестными и вычислены все значения искомых коэффициентов. Опытным путем удалось установить, что оптимальным является решение для N = 6 и и =5. При этом погрешность была минимальной:

-I <2х10-1°Дех>0. (11)

х

п /2

Г

х

с

с

с

2

X

п/2

Г

х

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

Следующим шагом было изучение гипергеометрической функции. Наиболее удобным выражением для моделирования является (4). Перепишем данную формулу для конечного числа членов ряда:

1FM(a,b,x) = £

/\

Л

J 1 J \AAaaJ

0 1,1 2,2 г,тт

Рис. 1. Радиальное распределение интенсивности ГГ моды с параметрами g=2, п=4 на расстоянии z=480 мм.

т %(Ъ)тт\' (12)

причем М должно удовлетворять условию:

^F1M (а, Ъ, х) - 1 Fl2M (а,Ъ, х)\ < е , (13)

где е определяется из (11).

Метод (12) с условиями (13) и (11) на первый взгляд удобный для компьютерного моделирования, оказался численно неустойчивыми при больших аргументах х (необходимыми для получения ГГ пучка). Имеется предел устойчивости вычислений комплексного распределения светового поля (6), при котором сумма (12) еще не расходится. Иллюстрация достижения такого предела приведена на рис. 1. Был проведен численный эксперимент, со следующими параметрами: выходное изображение имело 512x512 точек дискретизации, размер 4,4x4,4 мм, длина волны света I=633 нм, параметры моды g=2, п=4. Из рис. 1 видно, что на расстоянии г = 480 мм при достижении крайних значений аргумента происходит паразитный всплеск интенсивности. Это связанно с расходимостью

2,2 r,mm

Рис. 2. Радиальное распределение интенсивности ГГ моды с параметрами g=2, п=4 на расстоянии z=500 мм.

суммы (12) при больших значениях х (ниже будет показаны допустимые пределы).

Для (12) аргумент вычисляется по формуле:

ikr_ ' 2z

(14)

Таким образом, увеличение расстояния z уменьшает значение аргумента (14). И действительно, на рис.2 показана та же ГГ мода, но на расстоянии z=500 мм. В данном случае становятся заметны шумовые колебания интенсивности на краю, однако всплеска не наблюдается.

Используя (14) можно установить предельно допустимое значение аргумента функции (12). Так, используя самый точный тип вещественных чисел (long double) для a = 2,5-2i и b = 5 Re x < 48,02x103. Данный момент накладывает ограничение на использование выражения (4) при расчете гипергеометрических функций.

Рассмотрим выражение (3). При численном интегрировании оно примет вид:

1 F1(a, b, x) =

Г(Ь)

1

Г(а)Г(Ь - a) M

£ к M J I1 M

exp

xm M

(15)

причем М, аналогично как для (12) должно удовлетворять условию (13).

Исследования показали, что в этом случае не возникает ограничений, связанных с ростом аргумента гипергеометрической функции. Были получены хорошие результаты в тех случаях, когда выражение (12) неприменимо. На рис. 3 и 4 показана ГГ мода с теми же параметрами, что и на рис. 1 и 2, но на расстоянии г=300 мм.

0,968

0

1,1

b-a-1

a-1

2

2,01

1,34

0,67

Л

к ЛЛМллллл

1,1

2,2 г,тт

Рис. 3. Радиальное распределение интенсивности ГГ моды с параметрами д=2, п=4 на расстоянии z=300 мм.

Расчет аналитического решения (2) с помощью (12) и (15) сравнивался на основе среднеквадратичного отклонения:

5 =

1

К1 Е0( х, у, 2 )| -| Е (х, у, 2)|)

Е0( x, у, 2)|2

а) оШятмяшшеетж.

Рис. 4. Распределение инвертированной интенсивности (а) и фазы (б) ГГ моды с параметрами д=2, п=4 на расстоянии z = 300 мм

а)

0,484

(16) б)

2,2 г,тт

где Е0(х,у,2) - комплексная амплитуда поля, вычисленная с использованием (15); Е(х,у,2) -комплексная амплитуда поля, полученная при помощи (12).

В результате сравнения численной реализации (2) через сумму (12) и интегрирование (15), оказалось, что в области допустимых значений аргумента алгоритма (12) различие несущественно (точность до константы). На рис.5 продемонстрированы обе реализации на расстоянии z = 1000 мм. Из графиков видно, что решения идентичны и различаются только константой, поэтому для сравнения (12) и (15) с помощью (16) прибегали к нормировке обоих решений (максимальное значение приравнивалось к единице).

Результат моделирования распространения ограниченной апертурой Яе = 2,2 мм ГГ моды (7) с помощью интеграла (8) представлен на рис.6.

Рис. 5. Радиальное распределение интенсивности ГГ моды с параметрами д=2, п=4 на расстоянии z = 1000 мм:

а) реализация через интегральное представление (15);

б) реализация через полиномиальное представление (12)

2,2 г,тт

Рис. 6. Радиальное распределение интенсивности ограниченной апертурой ГГ моды с параметрами д=2, п=4 на расстояния z = 1000 мм

Из рис. 6 видно, что кривая интенсивности близка к кривой, изображенной на рис. 5, то есть с некоторой погрешностью вполне можно сформировать ГГ моду с помощью ограниченного апертурой распределения (7).

Оптимизация алгоритма

Вычисления гипергеометрических функций (полиномиальной или интегральной реализации) на компьютере для относительно больших значений параметра g является очень медленными. Например, для расчета моды с параметрами g=-10, п=4, при выходном изображении 512x512 точек дискретизации, раз-

0,694

0.208

0,347

0004

0

0,968

0

1,1

2 2 г,тт

0,074

0,037

0

1,1

мера 8x8 мм, длина волны света I = 633 нм, для расстояния г = 1000 мм на современном компьютере с частотой 2,2ГГц требуется 2 часа работы. Если рассчитывать моду при тех же условиях, но для g=-1, то потребуется всего 5 секунд. Для некоторых мод время вычисления занимало больше рабочего дня, поэтому желательно произвести оптимизацию алгоритма. Для этого был использован подход развертки радиально-симметричной функции, которая легко прослеживается в выражениях (2) и (6). Данный подход применялся в работе [19] и заключается в выявлении повторно вычисляемых близких по значению функций с различными значениями аргумента, сохранении и повторном использовании этих значений без затрат вычислительных ресурсов.

Перепишем (2) в следующем виде:

Еуп (г, р, г) = ^г, г, п, у) ехр(гпр) , (17) где для целого п

гу -1

^г,п,у) = ¿ИШ 2 х

ехр

1 Fl

1Л А I ■ 1ч

--(п\ - гу +1)

4

2 г

Г

п + гу+1 2

(18)

п +1 - гу , , , гкг

—-, п +1;-

2 11 2г

где зависит от знака индекса п:

(-1)'п, если п < 0 [1, иначе

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

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

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

Численные эксперименты

Решение уравнение (1) в цилиндрических координатах можно представить аналитическим решением (17), тогда мы получим комплексное распределение бесконечного светового поля и такие поля можно называть гипергеометрическими модами. Моделировать распространение светового поля (7) для конечной апертуры можно используя преобразование Кирхгофа. Далее рассматривается влияние на модовые свойства ГГ мод такого усечения светового поля.

Сравнение полученных распределений интенсивности и фазы на расстоянии г = 1000 мм при различных значениях параметра у приведено на рис. 8. Параметры расчета были выбраны следующие: изображения 2048x2048 точек дискретизации, размером 8x8 мм.

Для корректной оценки среднеквадратичного отклонения по формуле (16) полученные распределения интенсивности были нормированы по формуле:

Е '(х, у, г ) =

Е ( х, у, г ) тах|Е(х, у, г) .

(19)

Рис. 7. Метод радиальной развертки матрицы G

Выражение (19) приводит значения распределений поля по интенсивности в диапазон [0,1].

На рис. 9 показано сравнение распространения аналитической и ограниченной ГГ моды с параметрами у=-10, п=4, для входных картинок высокого разрешения (2048x2048), а на рис. 10 -зависимость среднеквадратичного отклонения от расстояния распространения, для более низкого (512x512), поскольку время вычисления возрастало с суток до недель.

Начальное падение СКО вызвано формированием моды, которая становится ближе к

п/2

Л

Параметры моды, СКО

Аналитическое решение

инверт. интенсивность

фаза

Моделирование с использованием интеграла (8)

инверт. интенсивность

фаза

г=-10, п=4

8 = 23,3%

у=0, п=4 8= 16,5%

^=10, п=4 8= 18,1%

О

Л\\

Рис. 8. Распределение интенсивности и фазы для ГГ мод на расстоянии вдоль оси распространения z = 1000 мм

г, мм 2600 5200 7800 10400 13000

Аналитическое решение © о О ¡ш Н

Моделирование с использованием интеграла (8) О о о О О

8, % 22,2 31,5 34,03 34,08 50,6

Рис. 9. Распределение инвертированной интенсивности для ГГ мод на различных расстояниях вдоль оси распространения z

аналитической, начиная с расстояния г « 2000 мм, что видно из рис. 10.

При формировании ограниченного апертурой ГГ-пучка его модовые свойства сохра-

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

5,%

50

40

30

2 600 5 200 7 800 10 400 13 000

Рис. 10. Зависимость среднеквадратичного отклонения от расстояния

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

гтах * , (20)

где Rх=Тхтг+Утг- размер входного Шо6-

ражения (для рис. 8 и 9 Rв = 4 мм).

Из (20) для ГГ моды, рассмотренной на рис. 9, получаем г « 5000 мм, что согласуется с результатами численных экспериментов. Из рис. 9 видно, что на г = 5200 мм пучок, полученный с использованием интеграла (8) имеет уже не такую четкую кольцевую структуру как аналитическое решение, т.е. можно сделать вывод что ГГ мода с этого расстояния начинает разрушаться.

Заключение

В заключении кратко сформулируем полученные результаты.

Для расчета аналитически полученного решения уравнения (1) в виде гипергеометрических мод (2) реализованы различные алгоритмы расчета.

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

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

Выполненено сравнение аналитического решения (1), которое получено для бесконечной ГГ моды, и моделирования распространения ограниченного поля (7) с помощью интеграла Кирхгофа.

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

Работа выполнена при финансовой поддержке российско-американской программы

"Фундаментальные исследования и высшее

образование" (грант CRDF RUXO-014-SA-

0б), а также гранта РФФИ 07-07-97600.

СПИСОК ЛИТЕРАТУРЫ

1. Siegman A.E. Lasers, University Science, Mill Valley, CA, 1986.

2. Miller W. Jr. Symmetry and Separation of Variables, Addison-Wesley Pub., MA, 1977.

3. Durnin J., Miceli J.J., Eberly J.H. Diffraction-free beams // Phys. Rev. 1987. V. 58.

4. BandresM.A., Gutierrez-Vega J. Ince-Gaussian beams/ // Opt. Lett. 2004. V. 29, No. 2.

5. Bandres M.A., Gutierrez-Vega J. Ince-Gaussian modes of the paraxial wave equation and stable resonators. // J. Opt. Soc. Am. A. 2004. V. 21. No.5.

6. Bandres M.A. Elegant Ince-Gaussian beams // Opt. Lett. 2004. V. 29. No.15.

7. Bandres M.A., Gutierrez-Vega J. Higherorder complex source for elegant Laguerre-Gaussian waves. // Opt. Lett. V. 29. No.19.

8. Schwarz U. T., BandresM.A., Gutierrez-Vega J. Observatiob of Ince-Gaussian modes in stable resonators // Opt. Lett. V.29. No.16.

9. Abramochkin E.G., Volostnikov V.G. Generalized Gaussian beams, J. Opt. A: Pure Appl. Opt., v.6, p.5157-5161 (2004).

10. Khonina S.N., Kotlyar V.V., ShinkaryevM. V, Soifer V.A., Uspleniev G.V. The phase rotor filter, J. Mod. Opt., v.39, no.5, p.1147-1154.

11. Beijersbergen M. W., Coerwinkel R.P. C., Kristiensen M., Woerdman J.P. Helical-wave front laser beams produced with a spiral phase plate // Opt. Commun. 1994. V.112.

12. Oemrawsingh S.S.R., Van Houwelingen J.A.M., ElielE.R., Woerdman J.P., Verstegen E.J.K., Kloosterboer J.G., Hooft G.W. Production and characterization of spiral phase plates for optical wavelengths. // Appl.Opt. 2004. V. 43. No.3.

13.SoskinM.S., Gorshkov U.N., VasnetsovM.V. Topological charge and angular momentum of light beams carrying optical vortices // Phys. Rev. A. 1997. V. 56. No. 5.

14. Котляр В.В., Хонина С.Н., Алмазов А.А., Сойфер В.А. Оптические чистые вихри и гипергеометрические моды. // Компьютерная оптика. 2005. № 27.

15.Котляр В.В., Скиданов Р.В., Хонина С.Н., Балалаев С.А. Гипергеометрические моды // Компьютерная оптика 2006. № 30.

16. Hart, J.F., et al. 1968, Computer Approximations/ New York: Wiley, 1968

17. Luke, Y.L. Mathematical Functions and Their Approximations. New York: Academic Press.

1975.

18. Hastings C. Approximations for Digital Computers (Princeton: Princeton University Press, 1955.

19. Балалаев С.А., Хонина С.Н. Реализация быстрого алгоритма преобразования Кирхгофа на примере бесселевых пучков // Компьютерная оптика. 2006. № 30.

CALCULATION OF HYPER-GEOMETRIC MODES

© 2007 S.A. Balalaev1, S.N. Khonina2, V.V. Kotlyar2

1 Samara State Aerospace University 2 Image Processing Systems Institute of Russian Academy of Sciences, Samara

An algorithm for calculating hyper-geometric modes of light field, which are the solution of a paraxial Schroedinger-type wave equation, is developed. Propagation of the hyper-geometric modes containing a phase spiral singularity (vortex) is computer simulated. Comparison is made between the rigorous analytical solution and an approximate solution derived using the integral Kirchhoff operator of propagation. Unlike the familiar paraxial (Gaussian and Bessel) modes, the major radius of the hyper-geometric modes increases with distance as a square root, with the distance between adjacent diffraction maxima (minima) being in inverse proportion to the radial coordinate.