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

РАСЧЕТ СИЛЫ ДАВЛЕНИЯ СВЕТА НА КРУГЛЫЙ ДИЭЛЕКТРИЧЕСКИИ ЦИЛИНДР С ИСПОЛЬЗОВАНИЕМ БЫСТРОГО ИТЕРАТИВНОГО АЛГОРИТМА И НА ОСНОВАНИИ АНАЛИТИЧЕСКОГО РЕШЕНИЯ

Котляр В.В., Налимов А.Г, Личманов М.А.

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

Аннотация

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

Введение

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

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

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

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

1. Расчет поля дифракции и силы давления света на диэлектрический цилиндр

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

у 2

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

1.1 Итеративный алгоритм

В [8] описан итеративный метод расчета дифракции произвольной электромагнитной волны, удовлетворяющей уравнению Гельмгольца, на цилиндрических диэлектрических микрообъектах с произвольной формой сечения (ТЕ-поляризация). Ядром итеративного алгоритма является уравнение Фредгольма второго рода:

Еи+1

где (2, у) е X ; Еп+1 и Еп - амплитуда электрического поля внутри области X на (п+1)-ом и п-ом шагах итераций; у - постоянная релаксации алгоритма, которая регулирует скорость его сходимости, _ 1к2 (е-1) 1 2п

р =-, к = — - волновое число света, е -

4 X

диэлектрическая проницаемость цилиндра, (§, п) е X - декартовы координаты внутри объекта,

Н0(1)(х) - функции Ханкеля 1-го рода нулевого порядка. В качестве начального приближения можно брать падающее поле:

Е,

Данный алгоритм сходится при условии кёе< 4п, где ё - диаметр сечения цилиндра. Для моделирования оптического захвата цилиндриче-

ских микрообъектов можно воспользоваться аналитическим решением дифракции гауссова пучка.

1.2. Аналитическое решение задачи дифракции

Рис. 1. Схема падения гауссова пучка с фокусом в точке (- 2о, Уо) на круглый цилиндр с центром в точке (0;0) Следуя [2], рассмотрим дифракцию двумерного непараксиального гауссова пучка на круглом однородном цилиндре (рис. 1). Для случая ТЕ-поляризации, когда (Ех, Ну, Н*) - отличны от нуля, напряженность электрического поля для непараксиального гауссова пучка можно записать в виде:

Ео х (г, ф) = Ео x (кг у

С =

Ю0л/п Г -—] ехр

к22д2

1 - д2 - гкдуо - т атс8т д

dд.

Напряженность электрического поля вне цилиндра равна сумме рассеянного и падающего полей:

Г

Ех = Ео X (кг У- +

Г

+Ео X Г«„СНЧкг) е'пф

где

а„ =

где к1 = —^ё^ - волновое число в среде с диэлек-

X

трической проницаемостью 81.

Связь между проекциями Н , Н* и Ех следует из уравнений Максвелла:

Н =-1 дЕк

у к д*

н2 = 1 дЕ±

к ду

1.3 Сила давления света на объект

Сила, действующая на объект со стороны падающего поля, рассчитанная в общем случае с помощью электромагнитного тензора Максвелла [9] в 2Б случае, может быть вычислена интегрированием по контуру £ вокруг объекта:

Ех = 0,

+Яе (НуН*) dy},

^ = 2ео|{2[IН* Г-81 Ех|2-|Ну

+Яе (Н*Н*) dz},

dz +

dy +

где 8о - диэлектрическая постоянная,

8о = 8.85.1о-12^-

о Н • м2

(в системе СИ).

2. Численное моделирование Исследуем работоспособность итеративного и аналитического методов расчета дифракции электромагнитной волны на микрообъектах. Далее, на основании результата вычислений поля дифракции электромагнитной волны на микрообъекте, рассчитаем силу давления света на этот микрообъект.

2.1. Сравнение итеративного и аналитического расчетов.

а) б)

Рис. 2. а) Амплитуда Ех поля дифракции непараксиального гауссова пучка на цилиндре, б) модуль проекции вектора Умова-Пойнтинга на ось распространения света 2.

Параметры: длина волны X =1 мкм, диаметр круглого цилиндра ¿1=1 мкм, диэлектрическая проницаемость цилиндра 8 =2, среды 81 =1, все поле 1,5x1,5 мкм, цилиндр помещен в центр перетяжки с радиусом ст =0,5 мкм

На рис. 2 представлено поле дифракции, рассчитанное по формуле (5). СКО между данным решением и расчетом картины дифракции при помощи алгоритма, описанного в [8], составляет о,о8% при параметрах: длина волны X = 1 мкм, диаметр круглого цилиндра ¿=1 мкм, диэлектрическая проницаемость цилиндра 8 =2, среды 81 =1, все поле 1,5х1,5 мкм, 256х256 отсчетов, объект помещен в центр перетяжки с радиусом ст =о,5 мкм, количество коэффициентов ряда Л/=8. Вектор Умова-Пойнтинга рассчитывается по формуле:

£ = 8оС2 [Е X Н ]

а его проекция на ось распространения света 2 будет равна:

£* =8оС2ЕхНу.

2

77 = — '-«

На рис. 3 приведены те же поля для цилиндра, смещенного на 0,25 мкм (вверх) по оси Y из перетяжки. По картине модуля проекции вектора Умова-Пойнтинга на ось распространения света хорошо заметно, что энергия в пучке отклоняется в сторону смещения цилиндра.

а) б)

Рис. 3. а) амплитуда Ex поля дифракции непараксиального гауссова пучка на цилиндре, б) модуль проекции вектора Умова-Пойнтинга на ось распространения света Z. Параметры те же, что и на рис. 2, цилиндр смещен на L=0,25 мкм по оси У

Точность расчета поля дифракции аналитическим методом зависит от количества коэффициентов ряда Сп (4). Зависимость невязки 5 между амплитудами электрических полей, полученных при помощи итеративного алгоритма и по формуле (5), от количества коэффициентов показана на рис. 4, рассчитываемая по формуле:

X (Ех (п, т) - Еп (п, т))2

Количество коэффициентее;

Рис. 4. Зависимость невязки 5 между амплитудами поля Ех, рассчитанного итеративным и аналитическим методами, от количества коэффициента Сп ряда

Количество коэффициентов ряда Сп изменяется

от -3..3 до -9..9, остальные параметры моделирования те же, что и для рис. 2.

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

тивным и аналитическим методами расчета дифракции от смещения круглого цилиндра вдоль вертикальной оси У (рис. 1) по линии, проходящей точно через перетяжку (7.=0).

а)

б)

Смешение 1_, мкм

Рис. 5. Зависимость невязки 5 между амплитудами

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

В аналитическом методе бралось: а) N=7 коэффициентов; б) N=9 коэффициентов

Как можно увидеть из рис. 5, невязка сильно возрастает при смещении цилиндра из фокуса в теневую область. Однако это можно исправить взятием большего количества коэффициентов Сп. Такая необходимость взятия большего количества коэффициентов Сп возникает из-за смещения значимых для расчета больших по модулю коэффициентов Сп от 0-го номера в сторону больших номеров, что проиллюстрировано на рис. 6.

На рис. 6а, б изображена зависимость модуля коэффициентов Сп от номера п. Как видно из графиков, при взятии того же количества коэффициентов и смещении цилиндра в вычислении принимают участие не все значимые коэффициенты, что дает увеличение погрешности расчетов.

Зависимость невязки 5 между картинами дифракции, рассчитанными итеративным и аналитическим методами, от смещения вдоль оси 2 дает иную картину. Так, на рис. 7 можно увидеть, что невязка незначительно возрастает при смещении цилиндра из фокуса вдоль оси распространения света на большее смещение, чем вдоль оси У. Параметры вычислительного эксперимента те же, что и для рис. 2.

Это объясняется другой зависимостью коэффициентов Сп от смещения, отображенной на рис. 8.

а)

б)

Рис. 6. Зависимость модуля коэффициента Сп от номера п: а) цилиндр совпадает с фокусом гауссова луча, б) цилиндр смещен из фокуса гауссова пучка на Ь=1 мкм по оси У

Рис. 7. Зависимость СКО между амплитудой Ех поля дифракции, рассчитанного итеративным и аналитическим методами, от смещения Ь цилиндра вдоль оси 2 через фокус при У=0

а коэффициента (

Рис. 8. Зависимость модуля коэффициентов Сп от

номера п. Перетяжка смещена от цилиндра по оси 2 на расстояние Ь=-1,56 мкм

Как видно из рис. 8, ширина полосы значимых коэффициентов Сп расширяется при смещении пе-

ретяжки из центра координат (и от центра цилиндра), и, чтобы это компенсировать, необходимо брать тем больше коэффициентов Cn, чем дальше смещена перетяжка гауссова пучка. В данном примере (рис. 8) взятие коэффициентов Cn от -9 до 9 дает невязку 5 между амплитудами Ex, рассчитанных аналитическим и итеративным методами, порядка 1% (рис. 7).

Взятие большего количества коэффициентов Cn сопряжено с увеличением вычислительной сложности и, следовательно, с увеличением временных затрат на расчет поля дифракции. Так, например, расчет поля дифракции размером 5x5 мкм (256x256 отсчетов) при параметрах: длина волны X = 1 мкм, диэлектрическая проницаемость среды е1=1, диэлектрическая проницаемость цилиндра е =2, диаметр цилиндра d=1 мкм (перетяжка расположена в центре координат) аналитическим методом при взятии коэффициентов Cn от -7 до 7 требует 40 секунд для компьютера с ЦП Celeron 1000, СКО между амплитудами Ex, рассчитанными аналитическим и итеративным методами составляет 0,038%. Расчет поля дифракции в том же случае, но при расположении перетяжки в координатах (Z=-1,56 мкм, Y=0 мкм) и взятии коэффициентов Cn от -13 до 13 требует 19 минут 53 секунды, невязка 5 между амплитудами Ex, рассчитанными аналитическим и итеративным методами составляет 0,077%.

2.2 Расчет силы давления

На рис. 9. приведен пример графиков зависимо -сти проекций силы, действующей на круглый цилиндр со стороны гауссова пучка, вдоль осей Y, Z от смещения цилиндра по тем же осям, из которых видна возможность «захвата» объекта лазерным излучением в оптическую ловушку.

а

Рис 9. Зависимость проекций силы, действующей на цилиндр, от его смещения вдоль той же оси: а) вдоль оси 2; б) вдоль оси У

Параметры моделирования: все поле размером 5x5 мкм, 256x256 отсчетов, длина волны X =1 мкм, диэлектрическая проницаемость среды 81=1, диэлектрическая проницаемость цилиндра 8 =1,3, диаметр цилиндра ¿=1,5 мкм, радиус перетяжки в фокусе ст =о,5 мкм. Линия смещения объекта на обоих графиках проходит через центр перетяжки. Из

рис. 9 видно, что сила, действующая на микроци-

и

линдр, равна по порядку величины 10-10 — при

м

мощности лазера гауссова пучка 100

мВт

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

На рис. 10 представлен график зависимости эффективности захвата Q от смещения вдоль осей У, 2, рассчитываемой по формуле:

Q = р

пР

где с - скорость света в вакууме, ^ - сила, действующая на объект, п =

те - показатель преломления объекта, Р - мощность излучения лазера.

а)

б)

Рис. 10. Графики зависимости эффективности захвата вдоль оси от смещения вдоль той же оси: а) ось распространения света 2; б) ось У

В случае силы, направленной вдоль рассматриваемой оси, Q>0, и Q<0 в обратном случае, что показывает возможность оптического «захвата» частицы. Физический смысл эффективности захвата -это часть излучения, пошедшего на создание силы, действующей на частицу.

Заключение В работе получены следующие результаты:

- проведено сравнение точности расчета дифракции гауссова пучка на диэлектрических микрообъектах аналитическим и итеративным методами

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

- оценена эффективность захвата диэлектрического цилиндра отдельно для смещения вдоль осей Y и Z (рис. 10) по формуле (10).

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

Работа поддержана российско-американской

программой «Фундаментальные исследования и

высшее образование» (BRHE), грант CRDF REC-

SA-014-02 и президентским грантом НШ-

1007.2003.01.

Литература

1. Zimmerman E., Dandliner R., Souli N. Scattering of an off-axis Gaussian beam by a dielectric cylinder compared with a rigorous electromagnetic approach. // J. Opt. Soc. Am. A, 1995, v. 12, p. 398-403.

2. Wu Z., Guo L.. Electromagnetic scattering from a multi-layerd cylinder arbitrarily located in a Gaussian beam, a new recursive algorithms. // Progress in electromagnetics research, 1998, PIER, v. 18, p. 317-333.

3. Ландау Л. Д. , Лифшиц Е.М. . Краткий курс теоретической физики. Механика. Электродинамика // Книга 1, М., Наука, 1969.

4. Lock J.A.. Calculation of the radiation trapping force for laser tweezers by use of generalized Lorenc-Mie theory I. Localized model description of an on-axis tightly focused laser beam with spherical aberration. // Appl. Opt., 2004, V. 43, p. 2532-2544.

5. Ganic D. , Gan X. , Gu M.. Exact radiation trapping force calculation based on vectorial diffraction theory. // Opt. Express, 2004, v. 12, no. 12, p. 2670-2675.

6. Nieminen T.A. , Heckenberg N.R. , Rubinstein-Dunlop H. . Computational modeling of optical tweezers. // Proceeings of SPIE, 2004, v. 5514, p. 514-523.

7. Pobre R. , Saloma C. . Radiation forces on nonlinear mi-crosphere by a tightly focused Gaussian beam. // Appl. Opt., 2002, v. 41, no. 36, p. 7694-7701.

8. Котляр В.В., Налимов А.Г., Скиданов Р.В. Быстрый метод расчета дифракции на цилиндрических диэлектрических микрообъектах // Компьютерная оптика, 2004, №25, с. 24-28.

9. Котляр В.В., Налимов А.Г., Скиданов Р.В. Расчет вектора Умова-Пойнтинга и силы давления электромагнитной волны на однородный диэлектрический цилиндр // Известия СНЦ РАН, 2004, т. 6, № 2.

10. Rohrbach A., Stelzer E.H.K. . Optical trapping of a dielectric particles in arbitrary fields. // J. Opr. Soc. Am. A, 2001, v. 18, p. 839-853.

м