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

УДК 535.42

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

© 2005 В.В. Котляр1, А.Г Налимов2, Р.В. Скиданов1

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

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

Введение

Обычно задачу манипуляции микрочастицами рассматривают для 3D случая, например, для частиц сферической формы. Впервые теоретическую задачу давления света на сферу с помощью теории Ми рассмотрел Де-бай в 1909 [1]. В общем случае сила давления света на микрочастицу должна рассматриваться с помощью максвелловского тензора напряжений электромагнитного поля [2]. Впервые сила давления слабосходящегося гауссова пучка на сферическую частицу с помощью теории Ми была оценена в [3], а для "острой" (непараксиальной) фокусировки гауссова пучка в [4,5]. Пользуются популярностью и приближенные методы расчета силы, действующей со стороны света на частицу: метод геометрической оптики с учетом фре-нелевских коэффициентов отражения и преломления [6], который применен для больших

значений параметра q = k0 а >> 1, где

С = ^ К 0 - л

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

ражение для силы давления света справедливо при условии q|n1 - п21 < 1, где п1 и п2 -

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

Для численного моделирования удобно перейти от 3D случая к двумерному случаю. Это приведет к существенному сокращению объема и времени вычислений, хотя задачи дифракции на сфере и на круглом цилиндре во многом аналогичны [11]. Поэтому в работе рассматривается дифракция цилиндрической электромагнитной волны ТМ-поляриза-ции на 2D однородном диэлектрическом цилиндре с сечением произвольного вида. Получены выражения для силы давления света на цилиндр. Вектор этой силы лежит в плоскости падения, перпендикулярной оси цилиндра, и выражается через векторы напряженности электрического и магнитного полей электромагнитной волны внутри и вне объекта.

Для расчета 2D поля дифракции ТМ-по-ляризованной волны разработан быстрый итеративный алгоритм, аналогичный алгоритму для ТЕ-поляризованной волны [12]. Алгоритм обладает работоспособностью при

I I к

условии, что ^п1 - п2| < 2 . Приводятся результаты численного моделирования дифрак-

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

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

В [1] приведена формула, выражающая сохранения полного импульса системы электромагнитного поля плюс объект V, ограниченный поверхностью S:

д$Ргё¥ + дР{)г ,

дг

дг

(1)

где Рг - координаты вектора импульса электромагнитного поля (V1 и S1 объем, ограничивающая его поверхность, который включает объект V еР1), который связан с вектором Умова-Пойнтига соотношением:

Р = 4=-!-

с 4%с

Е х Н

(2)

Р0г - координаты вектора импульса объекта,

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

2 _► 2

Е + Н

2 -5гк - ЕгЕк - НгНк ; (3)

дР

и1 ог

дг

света на объект;

1

4%

агк - максвелловский тензор напряжений

электромагнитного поля ( агк = аы ); Е, Н -векторы напряженностей электрического и магнитных полей в вакууме.

После усреднения по времени за период

т=— ш

монохроматического света:

Е(х, г) = ЯеЕ(х)вш },

Н(х, г) = Яе{Н(х)егшг} (4)

вместо уравнения (1) получим:

Р =/дРог

дг

= -§{агк)Пк^ , (5)

так как

д$ Р'п=0. (6)

\дг V IV

Можно показать, что:

/дР\

дг I 4%с

+ /яеЕ

— | Ке(гшЕу (Х)егш ' -гш

гшЕу (х)ег \Нг (х)ег + ^е(Еу (х)егшг )яе(гшНг (Х)егшг

^е(гшЕг (х)егшг )яе(Ну (х)ег ^е(Ег (х)егш )яе(гшНу (х)ег'

= 0, (7)

где Яе(...) - реальная часть комплексного

1 т

, (/(г)) = - $ I(г. Аналогично (7)

числа,

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

'Р , дг

'йр

к дг

= о

Для получения выражения для усредненного по времени тензора напряжений (3) учтем, что

(яе(Ег (х)егшг )яе(Е; (^х)егш')

1

= - Яе 2

; Е г (х)Е* (х)

Тогда вместо (5) получим:

(8)

Р = -1 $11

8% ' I 2

|2 I |2 I 12

Ех\ + \НХ\ -\Е \ -

- \Н„Г - |Ег|2 - |Н,|2^х +

у

+ )

Яе^Е* + НхН] у + Яе(ЕХ + НхН* )с152}

Р = - $11 у 8% $ 12

|2 I |2 , 12

Еу\ + Ну -|Ех| -

-|Нх|2-|Е\2 -|Нг|2 \пу + + Яе(ЕуЕ* + НуН* + Яе^Е* + НуН* К}

а гк =

^=-1 fí1

z 8л S [2

i2 i \2 \ 12 E\ + И -El -

2 2 2 - И I - IeJ - |и„| ysz +

y

y

+ .

Re(EzE; + И H* К + Re(EzE* + И H* К }

o*

F = i fji

У 8л Ц 2

\2 , i2 i 12 И,, - E - H

dSy +

+

Re(HyH* )dSz}

Fz = - fí1 Z 8лЦ2

2 2 2 Hz| -|E;| -И1

S +

И =- í OE;

y k Oz

H =

í OE;

k Oy

_ 2л

где к -— - волновое число света с длиной А

волны А. Аналогично (10) сила давления света с ТМ-поляризацией на 2D объект будет иметь следующие проекции

(Ех -Ну -Нг - 0):

(9)

O Oz

где dS; = -^dxdy , dSy = ~yd;dy ,

F = -1 fJ 1 y 8л f I 2

F = 0

2 2 2 E - E - И

dSy +

dSz - ^х^у, Е1 - Ех, Е2 - Еу, Е3 - Ег (для

Hi и Ег. аналогично).

Перепишем выражения (9) для силы давления света на микрообъект в 2D случае. Для

ТЕ-поляризации ( Нх - Еу - Ег - 0 ) электрическое поле направлено вдоль оси X: Ех ^ 0, Z - оптическая ось, 2D-объект имеет вид цилиндра с произвольной формой сечения и имеет бесконечную протяженность вдоль оси X. Плоскость YOZ - плоскость падения света. В этом случае соотношение (9) примет вид:

К - 0,

+

ReEyE* К}

F = ^ fí 1

8л-! I 2

2 2 2 EÁ - Ey - |И;|

dS +

(12)

+ Re{HzH*y ^у }. (10)

Здесь £1 уже контур, охватывающий сечение

объекта в плоскости YOZ. Сила - направлена вдоль оптической оси и является аналогом рассеивающей силы для рэлеевских частиц [7], а Ку - направлена поперек оптической оси и является аналогом градиентной силы [7].

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

+ RefeEy )dSy }, где (как и в уравнении (10))

dSy = nydl = sin tydl = dz

и

dSz = nzdl = cos <$dl = dy. dl - элемент дуги.

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

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

задачу дифракции [13] (ф - Нг).

Ф( ^ z) = Фо( ^ z)s-

íko2(s-1)

4

■ff<KS, л) И

Ч)2 + (z-л)2 )d;dn +

+ ^f л)И02)(kV(y Ч)2 + (z-л)2 )й

(13)

внутри объекта площадью V, ограниченного контуром S, т.е. (y,z)eV и

(11)

ф( У, z) = фо( У, z) -

1)

4s

■{{ф(£, л) Н

Ч)2 + (г -л)2 ^л

+

+ ^ | ф(£, л) Н 02) (у -£)2 + (г-л)2 У

дп

(14)

снаружи объекта при (у, г) £ V, где в - диэлектрическая проницаемость внутри объекта, Н 02) (х) - функции Ханкеля нулевого порядка,

ф 0 - амплитуда падающей световой волны, П - внешняя нормаль к поверхности объекта.

Интеграл по контуру в (13) и (14) можно представить через преобразование Фурье, введя функцию контура объекта.

'1, (У, г)е 5

G(У, г) = ■

0, (У, г)£ 5 •

(15)

Тогда

д

,(У, г) — Н 02) \Ч (у -£)2 + (г-л)2 Щ =

:ЛЗ-1 [з(ф( у, г )в( у, г)Пу )■

•3

2 2 2 +л2

Н 12)

22 +л

+

+ 3(ф( у, г )G( у, г )пг )■

■3

^^ Н1(2) ([^Ты)

л/£2 +л2 ,

(16)

(пу , пг )

где п = \пу, п2) - единичный вектор внешней

■3

(Ф „ (£, л)р(£, л))■ з(н02) (^^Т^^тл

-у 3-1 [з(ф „ (£, л^ (£, л)Пу )■

■3

.V«

2 , „2 +л

Н

(2)

22 2 + л2

3-1 [з(ф „ (£, лМ£, л)Пг )■

■3

N

.лй2 +л2

внутри объекта (у, г) е V и

фп+1( У, г) = уф0( ^ г) -у

+(1 -у)ф п(У, г) (18)

?^02(£-1) ^ 4в

■3

3(ф„ (£, л) р(£, л)>3(н 02) (кЛ

2 , „2 +л

3-1 [3(ф „ (£, лМ£, л)пу )■

■3

Н 12) (кдй2^)

3-1 [3(Ф „ (£, лМ£, л)пг )■

4в

• 3

Я

2 2 2 +'

+ (1 - 7)Фп (у, г)

(19)

нормали к контуру S, (у, г) е 5, 3 и 3-1 -

прямое и обратное 2D Фурье-преобразования. г / \ „ , ,

. снаружи объекта (у, г)£V , где ф, и ф -

Аналогично, интеграл типа свертки по площади объекта V в (13) и (14) можно перепи- рассчитанные комплексные амплитуды поля сать через преобразование Фурье с помощью функции объекта

, ч Г1, (£, л)е V

на (п+1)-й и (п)-й итерациях, у - параметр релаксации.

Численное моделирование

Для проверки работоспособности алго-

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

Итеративный алгоритм решения интегрального уравнения (13) будет иметь вид:

ф п+1( У, г) = уф0( У, г)в-у

л о (в -1) 4

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

£

£

£

30 -

25

20 ■■■ 5Й ' О 15

(J

10 5 0 -

-

5 10 15 20 25 30 35 40 45 50 55 60 65 70 Ы итерации

5 10 15 20 25 30 35 40 45 50 55 60 65 70 Ы итерации

б)

Рис. 1. а - среднеквадратичное отклонение модуля амплитуды магнитного поля, рассчитанного данным алгоритмом, от точного решения в зависимости от числа итераций; б - среднеквадратичное отклонение двух модулей амплитуд магнитного поля, рассчитанных данным алгоритмом, для двух соседних итераций, в зависимости от номера итерации

Если сравнить данные графики (рис. 1) с аналогичными, полученными для ТЕ поляризации [12], можно заметить, что устойчивость алгоритма для ТМ поляризации выше, чем для ТЕ, но точность ниже. Это обусловлено ошибкой расчета интеграла по контуру объекта на границе раздела двух сред.

Параметр релаксации у был подобран экспериментально для случая дифракции на круглом цилиндре. Подбор был осуществлен на поле 256x256 отсчетов, 2,5x2,5 мкм. Диаметр цилиндра был равен длине волны в d = X=1 мкм, диэлектрическая постоянная цилиндра, находящегося в вакууме, в =2. При этих параметрах среднеквадратичное отклонение рассчитанного за 13 итераций модуля амплитуды напряженности магнитного поля от точного решения в зависимости от параметра у показано на рис. 2.

Как видно из рис. 2, достаточно широкий диапазон параметра у приводит к устойчивой работе алгоритма, наилучшим являет-

S 0 3 0. 35 0 4 0. 45 0 5 0. 55 /

Рис. 2. Зависимость невязки на 13-ой итерации от коэффициента у

ся у =0,35..0,4. Экспериментально было найдено, что алгоритм не расходится при усло-

вии, что qn1

I к 2к

- ni\< 2, где q = Y a-

a - ра-

диус окружности, в которой вписано сечение цилиндрического объекта, n2 и n1 - показатели преломления объекта и среды.

На рис. 3 приведены примеры расчета картин дифракции на круглом цилиндрическом диэлектрическом объекте непараксиального гауссова пучка с ТЕ- и ТМ-поляризован-ной волной. Параметры расчета: все поле 10x10 мкм, 512x512 отсчетов, длина волны X=1 мкм, диаметр цилиндра D=1 мкм, диэлектрическая проницаемость объекта s =2, диэлектрическая проницаемость среды s 0=1, радиус перетяжки пучка в фокусе а =0,5 мкм. На рис. 3а,б объект расположен точно в фокусе, на рис. 3в,г смещен вверх (по оси Y) на половину длины волны. Время расчета поля с ТЕ поляризацией на компьютере PIII Celeron 1000 MHz составляет 31 секунду, всего 155 итераций. Время расчета поля с ТМ поляризацией составило 1 минуту и 3 секунды, было сделано 72 итерации. Критерий остановки итераций в обоих случаях одинаковый: невязка между шагами должна быть меньше 0,0001.

Вектор Умова-Пойнтинга

1

S =

4кс

Re

E х H

Рис. 3. Модуль напряженности магнитного поля (а,в) ТМ-волны и электрического поля (б,г) ТЕ волны, которые описывают дифракцию непараксиального гауссова пучка на круглом диэлектрическом цилиндре: цилиндр в фокусе пучка (а,б) и смещен вверх на пол длины волны (в,г)

для 2D задачи имеет вид: для ТЕ-поляризации:

5 = -

1

4л:ск

еу 1т

К,

dE*

_х

dу

(

+ е, 1т

К,

dE

х

dz

, (20)

где еу и е - единичные векторы вдоль осей

Y и Z, 1т(...) - мнимая часть числа, с - скорость света в вакууме; для ТМ-поляризации:

5 =

1

4л:с

(

е 1т

у

Н

dу

А

(

+ е 1т

Н

dH

V

dz

, (21)

На рис. 4 показаны поля направлений вектора Умова-Поинтинга (20) и (21), наложенные на картины дифракции плоской (рис. 4 а,б) волны, и гауссова пучка (рис. 4 в,г) с длиной волны X=1 мкм, диаметр цилиндра D=1 мкм, диэлектрическая проницаемость объекта £ =2, диэлектрическая проницаемость среды £ 0=1, радиус перетяжки пучка в

фокусе а =0,5 мкм.

Амплитуда падающего непараксиального пучка рассчитывалась по формуле [14]:

Рис. 4. Поле направлений вектора Умова-Пойнтинга для дифракции электромагнитной ТЕ-поляризованной (а,в) и ТМ-поляризованной (б,г) волны на круглом диэлектрическом цилиндре: освещающая волна плоская (а,б), гауссовая (в,г) (цилиндр в перетяжке гауссова пучка)

EZ 0( ^ у) =

1

а 0л/2л

I ехр|" Т^Г + Й + Ул/1 2 , (22)

-ж [ 2а 0 ]

где а 0 - параметр гауссова пучка, связанный

с радиусом его перетяжки ю соотношением: =

кю '

Силы, расчитанные по формуле (10), действующие на круглый цилиндр в случае ТЕ-поляризованной волны, представлены на рис. 5. Параметры эксперимента: все поле 10x10 мкм, 512x512 отсчетов, длина волны X = 1 мкм, диаметр цилиндра D=1 мкм, диэлектрическая проницаемость объекта £ =2, диэлектрическая проницаемость среды £ 0=1, радиус перетяжки пучка в фокусе а =0,5 мкм. Смещение от центра L имеет размерность мкм. Силы даны в безразмерных единицах.

Аналогичные силы в случае ТМ-поляри-зации, рассчитанные по формуле (12), представлены на рис. 6.

L

а)

0,15 0,1

0,05

0

-0,05

-0,1 -0,15

б)

Рис. 5. ТЕ-поляризация: а - зависимость силы Fz от смещения L объекта вдоль оси Z через центр перетяжки ^=0); б - зависимость сил Fz и FУ от смещения L объекта вдоль оси Y через центр перетяжки ^=0)

Заключение

В работе получены следующие результаты: разработан итеративный алгоритм для быстрого расчета поля дифракции ТМ-по-ляризованной электромагнитной волны на цилиндрических диэлектрических микрообъектах;

с помощью быстрого преобразования Фурье за 74 итераций, которые на компьютере Celeron 1000 MHz выполняются за 1 минуту 18 секунд, было рассчитано поле дифракции плоской волны на круглом цилиндре, отличающееся от точного решения в среднем на 4%;

численно показано, что параметр релаксации у итеративного алгоритма можно варьировать в широком диапазоне 0,3< у <0,525 для достижения практически одинаковой точности решения задачи дифракции;

рассчитано поле дифракции непараксиального гауссова пучка на диэлектрическом круглом цилиндре, диаметр которого сравним с диаметром перетяжки и расположен вблизи этой перетяжки; рассчитан вектор Умова-Пойнтинга для дифракционной картины ТМ- и ТЕ-поля-ризованной волны на диэлектрическом

0,07 0,06 0,05

ш

г 0,04

Ё о,оз 0,02 0,01 о

а)

-0,05

б)

Рис. 6. ТМ-поляризация: а - зависимость силы Fz от смещения L объекта вдоль оси Z через центр перетяжки (Y=0); б - зависимость сил F и F от смещения L объекта вдоль оси Y через центр перетяжки ^=0)

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

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

Работа выполнена в рамках российско-американской программы "Фундаментальные исследования и высшее образование" (грант CRDF REC-SA-014-02), поддержана грантом НШ-1007.2003.1, грантом Министерства образования и Правительства Самарской области.

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

1. Debye P. Der Lichtdruck and Kugeln von beliebige Material // Ann. Phys. 1909. V. 30.

2. Ландау Л.Д., Лифшиц Е.М. Механика. Электродинамика. М.: Наука, 1969.

3. Kim J.S., Lee S.S. Scattering of laser beams and the optical potential well for a homogeneous sphere // J. Opt. Soc. An. 1983 V. 73.

4. Barton J.P., Alexander D.R., SchaubS.A. Theoretical determination of net radiation force and torque for a spherical particle illuminated by a purged laser beam // J. Apl.

Phys. 1989. V. 66.

5. Ren F., Grehan G., Gonesbet G. Radiation pussure forces exerted on a particle located arbitrarily on a Gaussian beam by using the generalized Lorents-Mie theory, and associated resonance effects // Opt. Commun. 1994. V. 108.

6. Ashkin A. Forces of a single-beam gradient laser trap on a dielectric sphere in the ray-optics region // Biophys. J. 1992. V. 61.

7. Navade Y., Asakure T. Radiation forces on a dielectric sphere in the Rayleigh scatterry regime // Opt. Commun. 1996. V. 124.

8. Lemire T. Coupled-multipole formulation for the threatment of electromagnetic scattery by a small dielectric particles of arbitrary sphere. //J.Opt.Soc.Am. A. 1997. V. 14.

9. Rohrbach A., Stelzer E.H. Optical trapping of dielectric particles in arbitrary fields // J. Opt. Soc. Am. A . 2001. V. 18.

10. Pobre R., Salome C. Radiation force on a nonlinear microsphere by a lightly focused Gaussian beam // Appl. Opt. 2002. V 41, №» 36.

11. Kozaki S. Scattering of a Gaussian beam by a homogeneous dielectric cylinder // J. Appl. Phys. 1982. V. 53.

12. Котляр В.В., Налимов А.Г., СкидановР.В. Быстрый алгоритм расчета дифракции электромагнитной волны на цилиндрических диэлектрических объектах // Компьютерная оптика. Самара: ИСОИ РАН. 2003. Вып. 25.

13.Kotlyar V.V, Lichmanov M.A. Analysis of light diffraction by micro-optics using finite elements method // Opt. Mem. and Neut. Net. 2001. V. 10. No. 4.

14. Petersson L.E., Smith G.S. Three-dimensional electromagnetic diffraction of a Gaussian beam by a perfectly conducting half-plane // J. Opt. Soc. Am. A. 2002. V.19. No.11.

CALCULATING OF THE UMOV-POITING VECTOR AND PRESSURE FORCE

OF THE ELECTROMAGNETIC WAVE ON THE HOMOGENEOUS DIELECTRIC CYLINDER

© 2005 V.V. Kotlyar1, A G. Nalimov2, R.V. Skidanov1

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

Algorithm for calculating of diffraction of the electromagnetic monochromatic wave on the dielectric micro objects with the sizes same to wavelength is considered. Comparisons of the numerical modeling of diffraction of TE- and TM- polarized waves on a micro objects with exact solving of this problem are presented. Using this algorithm Umov-Pointing vector and pressure force of the non-paraxial cylindrical Gaussian wave on the circle micro cylinder are calculated.