CC BY
40
Вычислительные технологии Том 10, Специальный выпуск, 2005
АНАЛИТИЧЕСКИЕ АЛГОРИТМЫ УСРЕДНЕНИЯ В ЗАДАЧАХ ДИФРАКЦИИ СВЕТА НЕСФЕРИЧЕСКИМИ ЧАСТИЦАМИ*
Л.Е. Парамонов, В. А. Шмидт, Г. В. Черкасова Красноярский государственный технический университет, Россия e-mail: param@ktk.ru, morgen@akadem.ru
T-matrix method and representation theory of the rotation groups are applied to develop analytical methods for computation of the optical characteristics of ensembles of the nonspherical particles, namely elements of the scattering matrix and scattering flux for any solid angles.
Введение
Знание и моделирование оптических характеристик необходимы во многих областях науки и техники: для интерпретации данных дистанционного зондирования атмосферы и океана [1, 2], идентификации биологических клеток [3], математического обеспечения приборов оптической диагностики — спектрофотометров [4], нефелометров поляриметров Сток-са [5], проточных цитометров [3].
1. Вектор Стокса и матрица Мюллера
Пусть направление распространения излучения задается вектором п = (в, ф), где в, ф — углы сферической системы координат. Вектор напряженности электрического поля имеет вид
Е = Ев ев + Е^ е^, (1)
где е^, е^ — орты сферической системы координат.
Представление оптических характеристик в базисе [6]
е-1 = —=(ee - e+i = —f(e® + ieV (2)
будем называть СР-представлением, а орты е-1, е+1 интерпретировать как право- и лево-циркулярно поляризованные излучения единичной интенсивности.
* Работа выполнена при частичной поддержке Российского фонда фундаментальных исследований (грант № 04-05-64390).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.
В дальней волновой зоне (кг ^ 1) рассеянная частицей волна сферическая [7]:
Е-1 Е+1
егкт
Е-1 Е+1
(3)
где С = {Орд}рл=+1—1 — амплитудная матрица рассеяния в СР-представлении; г — расстояние до точки наблюдения; к — волновое число; индексы в, г означают рассеянное и падающее поля соответственно.
Определим параметры Стокса в СР-представлении следующим образом:
1-2 = Е-1Е+1 = 2(3 - ги), 1-о = Е-1Е-1 = 1(/ + V),
/+о = Е+1Е+1 = 2(1 - V), 1+2 = Е+1Е-1 = 2(3 + ги). (4)
Здесь /, и, V — параметры вектора Стокса в ЬР-представлении [8].
Такой выбор параметров Стокса в СР-представлении отличается от традиционно используемых [9] и обусловлен необходимостью согласования параметров Стокса и элементов матрицы Мюллера — преобразование параметров Стокса описывается матрицей Мюллера и имеет в СР-представлении следующий вид:
/р-рЫ = 0рд ОЬ1гд-д(пг),Р,д,Р,д= -1, +1, (5)
где знак нижнего индекса /р-р совпадает со знаком р.
Использование теории представлений групп [10] и канонических базисов в пространстве решений векторного уравнения Гельмгольца и пространстве И3 позволяет разделить переменные в выражениях для элементов амплитудной матрицы рассеяния в СР-представлении в терминах элементов Т-матрицы [11]:
1 тете п п'
Орд =2ЕЕЕ Е 1пп' (-1Г+™' Атт' (^гКт^К!' (6,)^
п=1 п'=1 т=-пт'=-п'
Т(рд) = Т11 + пТ 12 + РТ 21 + РпТ 22 (6)
тпт'п' тпт'п' 1 У тпт'п' 1 I тпт'п' 1 гЧ -1 тпт' п'' "ч/
Атт' (<Рэ,&) =
Ьап' = гп'-п-1[(2п + 1)(2п' + 1)]1/2,р, п = -1,1.
Здесь ¿птт,(6) — функции Вигнера [6]; (6г,^г), (68,ф8) — направление соответственно падающей и рассеянной волн.
Факторизация по параметрам падающего, рассеянного излучений и ориентации частиц является основой для разработки эффективных аналитических алгоритмов решения классов задач, связанных с оценкой: 1) потоков рассеянного излучения в произвольных конических телесных углах (усреднение интенсивности рассеянного излучения); 2) оптических характеристик ансамблей несферических частиц произвольной ориентационной структуры (ориентационное усреднение); 3) при некогерентном падающем излучении произвольной структуры (направление, интенсивность и поляризация) (усреднение по параметрам падающего излучения) [11, 12].
2. Метод Т-матриц
В настоящее время метод Т-матриц [13] — один из самых эффективных методов решения задач дифракции электромагнитного излучения частицами несферической формы. Если следовать формализму метода Т-матриц, падающее и рассеянное поля представляются соответственно [13, 14] как
те п
ЕЧГ) ^Е Кп:^Мтп(кг) + 6тпК^тп(кг)], (7)
п=1 т=—п те п
Е'(г) = ^ ^ [ртпМтп(кг) + Ятп^тп(&г)].
п=1 т=—п
Здесь Мтп(кг), Ктоп(кг), RgMmn(kг), RgNmn(kг), п = 1, 2,..., — п < т < п, — канонический базис в пространстве решений векторного уравнения Гельмгольца [15], преобразующийся по неприводимым представлениям веса п группы вращений и не зависящий от способа задания вращения.
Коэффициенты разложения связаны линейным преобразованием [13]
Р _ Ч
так называемой Т-матрицей — инвариантом относительно направления распространения падающего излучения в фиксированной системе координат. Т-матрица зависит от размера, формы, внутренней структуры и ориентации частицы, а также от ее относительного показателя преломления.
Пусть g — вращение системы координат (Р ^ Ь), тогда Т-матрица преобразуется согласно [15] следующим образом:
п п'
Тпгпт'п'(Ь) = Е^ Е^ Пттп1 ^ )Тт^тт2п' (Р)^етот)2(g ), ^ = 1, 2, (9)
т\ = — п т2 = —п'
или в матричной форме Т(Ь) = D(g—1 )Т(Р-1), где g — последовательное вращение системы координат относительно подвижных или неподвижных осей; П^п^ (g) — соответствующие элементы унитарной матрицы неприводимого представления веса п группы вращений; Ь, Р обозначают соответственно лабораторную и связанную с частицей физическую системы координат. В дальнейшем используется последовательное вращение системы координат относительно подвижных осей Z, У, Z, для которого 1) = П^т (авт) —
П-функции Вигнера [6], а, в, 7 — углы Эйлера. При вращении пространства в формуле (9) —1
g 1 следует заменить на g.
Полагаем, что Т-матрица частицы в физической системе координат известна. В дальнейшем ограничиваемся рассмотрением осесимметричных частиц, ось вращения которых совпадает с осью Z физической системы координат.
3. Матрица рассеяния монодисперсного ансамбля несферических частиц
Введем лабораторную систему координат (Ь) таким образом, чтобы направление падающей волны совпадало с осью Z, а плоскость рассеяния — с плоскостью XZ, что соот-
Т11 Т12 а
Т21 Т22 Ь
(8)
Рис. 1. Лабораторная (Ь) и физическая (Р) системы координат.
ветствует 6г = ^г = = 0. Углы Эйлера а, в, 7 определяют последовательное вращение системы координат (Ь) к системе координат (Р) относительно подвижных осей Z, У, ^ (рис. 1).
В этом случае матрица Мюллера является матрицей рассеяния, а элементы амплитудной матрицы будут равны
тете п
Орд = -2ЕЕ Е гп'-п[(2п + ^ + ^—Т+'^б*)^(Ь),р,п = ±1. (10)
п=1 п'=1 т=-п
Для монодисперсных ансамблей частиц с заданным распределением по ориентациям р(а, в, 7) элементы матрицы рассеяния в лабораторной системе координат будут иметь вид
л2п лп л2п
(ОрдО*т) = йа I в1пвdв ¿^р(ав^)Орд(68,ав1)О*т(68,ав7),Р,П,Р,П = ±1. (11) ./0 Jo
Численное усреднение элементов матрицы рассеяния не всегда эффективно, поскольку требует больших вычислительных затрат. В ряде случаев эффективен подход с использованием разложения элементов матрицы в ряды по обобщенным сферическим функциям [7, 10] или функциям Вигнера [6].
Коэффициенты разложения являются компактным и удобным способом хранения информации. При известных коэффициентах расчет элементов матрицы рассеяния может быть произведен с минимальными вычислительными затратами, а эти данные могут быть многократно использованы при решении задач однократного рассеяния, а также уравнения переноса поляризованного излучения в атмосферах планет [9].
3.1. Хаотически ориентированные частицы р(а@7) =
1
8п2
В [16] получены аналитические формулы для коэффициентов разложения элементов матрицы рассеяния в ряды по функциям Вигнера для хаотически ориентированных частиц:
(ОрдОрд) 'У ] Ур-р)д-д"'р-р)д-д\^ 5
п1=таж(|р-р|,|д-д|)
9%- р.д-дС- г,.д-а(65),Р,П,р,П= -1, +1.
(12)
те
Коэффициенты разложения р й(^) в формуле (12) имеют вид
р—р,^
п+п! / 2 Л + 1 \ 1/2 тгп(п , п+д-д)
2П + 1 \ с пр с пт п(РЯ'5'')
р——У^'"1 1 / у / у \ 2п + 1 у спрп1р—р / у с»ггп,п!д—<гПтп?тг»г'
«п! = (2П + 1) 2 ' + М Сп» С пт П
Ур—р,я—я У^"-1 ^ V / у / У I 2п + 1 / прп!»—р / у иптп1?—я '
п=1 п=таж(1,|п—п!|) ^ т=—тгп(п,п+д—д)
(13)
т = т + ^ — д, р, д, р, ^ = —1, +1,
где С'п71тп1т' — коэффициенты Клебша — Гордона [6].
В частности, для частиц, обладающих осевой симметрией, запишем уравнения
те
П(ряря) = у^ (2 +1)в(р«) в(ря)* (14)
тптп / у \ I"1 1 / тппц тпп1' V /
п! = |т—
п+п1
в (ря) = сп'9 А(р?)
втпп1 / у сптп1д—тАпп'п1'
п'=таж(1,|п—п! |)
'
п —п
тт(п,п')
а(ря) = ^__с п'т! Т
пп'щ = 2(2п' + У Спт1п10 г
для сферических частиц
(р?)
т1 пт1 п
т! =—тт(п,п')
птрпртп = ^¿т^В«®*, Вря) = ^п + 1)1/2(«п + мЬп), (15)
где ап, Ьп — коэффициенты Ми [17]. Отметим ряд важных свойств симметрии, используемых при численной реализации алгоритма в случае осесимметричных частиц:
Т(ря) = х т(Р?) Т(Р?) = Т(—Р—А(Р?) = (_1)п+п' +п! А(—р—я) в(ря) = в(—р—я)
Т тпт'п' хтт' Ттпп', Т тпп' Т — тпп' ' лпп'п1 ( 1) лпп'п1 , втпп1 в—тпп!,
тт(п,п')
А(ря) =
п! (2п' + 1)1/2
1Т 1(ря) С п0 + \ л с пт! Т к(ря)
2 Т 0пп' сп'0п10 + / у сп'т1п10Т т1пп'
(16)
т!=1
к =1, если п + п' + п1 — четное число, иначе к = 2,
Т 1(р«) = (т 11 + РУТ22 ) Т2(Р?) = (УТ12 + РТ21 ) = -1 +1
Тт1 пп' \Тт1пт1п^ руТт1пт1п'/, Т т!пп' \уТт1пт1 п' "г" рТт1пт1п'/, р У , Т-1-.
Полученные формулы (13)-(16) имеют компактный вид, удобный при численной реализации. Для осесимметричных частиц алгоритм был численно реализован, в частности, для сфероидальных частиц получено совпадение всех значащих цифр с результатами, представленными в работе [18].
3.2. Аксиально ориентированные частицы 7) = 4—^Р(в)
Под аксиально ориентированными частицами будем понимать ансамбли осесимметричных частиц, ориентационная структура которых в лабораторной системе координат не зависит от углов а, 7. Вид функции р(в) зависит от размера, формы, удельного веса вещества
частиц, вязкости окружающей среды и физически может быть реализован под действием ориентирующих факторов — геофизических полей (гравитационных, электрических, магнитных, поля скоростей и др.).
В этом случае в формулах (13), (14) величина имеет вид
ОрФФ = /Ырд)в(м)*\ (17)
те тгп(п,п')
В»в) = 2 Ё + !)1/2 Ё ^(в)ТЙпп'(в),
п'=1 т,1 =—тгп(п,п')
скобки означают ориентационное усреднение. Отметим свойство симметрии, используемое при численной реализации:
В--^ = (-1)™-« ВЫ. (18)
В случае р(а,в, 7) = —-$(еов в — сов во) в формуле (17) отсутствует ориентационное
4п2
усреднение по в.
4. Потоки рассеянного излучения в произвольных конических телесных углах
Пусть рассеянное частицей излучение попадает на приемник, характеризуемый в лабораторной системе координат углами ($, ф) и линейным углом конического телесного угла 2$0. Определим поток рассеянного излучения, проходящего через поверхность приемника.
Введем лабораторную систему координат (Ь) таким образом, чтобы направление падающей волны совпадало с осью Z, и физическую систему координат (Р), связанную с частицей (рис. 2). Углы Эйлера (а, в, 7) определяют ориентацию частицы и, соответственно, поворот системы координат (Ь) ^ (Р).
Пусть система координат (Б) связана с приемником, а поворот системы координат (Ь) ^ (Б) относительно подвижных осей тогда определяется углами Эйлера (ф, 0).
Рис. 2. Лабораторная (Ь), физическая (Р) системы координат и система координат (5), связанная с приемником.
Соотношения T-матриц в системах координат (L), (P), (S) имеют вид
T(L) = Б(ав7 )T(P )D-1 (а^т),
T(L) = D(0^O)T(S)D-1(0^O), (19)
T(S) = D-1(^0)D(a,5Y) • T(P) ■ D-1(a,5Y)D^0). Поток излучения через площадь приемника в системе координат (S) равен [19]
Ф = / (I+o + I- o)r2du = (/+0) + </- 0), (20)
п
где du = sin 0sd0sd^s; О = ((0s,^s)| 0 < ^s < 2п, 0 < 9s < $0), (9s,ps) — направление рассеянного излучения в системе координат (S).
Так как элементы вектора Стокса I+0,1-0 инвариантны относительно вращения плоскости референции, нахождение потока излучения сводится к усреднению элементов матрицы рассеяния:
<I±0> = £ <C±1q C±1q)/J-,. (21)
q,<j=-1,+1
Здесь Iq-q — элементы вектора Стокса падающего излучения относительно плоскости ^ = ф в лабораторной системе координат.
Поскольку в системе координат (S) направление падающего излучения определяется сферическими углами ($, — п), элементы амплитудной матрицы можно записать в виде
..те n
Cpq(6s,Vs) = (—i)meim^dnm(^b^W,
n=1 m=-n
те
B[mn?($) = Y. tnn' 4m¡n> w, (22)
n'=1
^тРПпПп' W = E ^m' WТтПт'п' (S)'
m'=-n'
Отметим, что элементы Cpq зависят также от параметров •
После аналитического усреднения в (21) выражения элементов матрицы рассеяния будут иметь вид
n+n
(cPqc;4) = ¡J]Bmq)mBmS*m Y. (—lY^cnpn-pCnmn-midno ш>, (23)
пПт n1 = |n-ri|
где
„,&0 ( 1 — cos $0, n1 = 0,
(dno (0s)> = dno (^s)sin ^s = ] sin dnii (^0), ni > 0.
0 I Vni(ni + 1) 01 '
Оценка потоков рассеянного излучения изотропными ансамблями частиц рассматривается в [12], в том числе при некогерентном падающем излучении различной геометрии [20].
п
Заключение
Аналитические алгоритмы могут быть использованы для разработки математического
обеспечения приборов оптической диагностики, а расчеты на их основе — для интерпретации экспериментальных данных.
Список литературы
[1] Lee Z.P., Carder R.L., Marra J. et al. Estimation primary production of depth from remote sensing // Appl. Opt. 1996. Vol. 35. P. 463-474.
[2] Зуев В.Е., Кабанов М.В. Оптика атмосферного аэрозоля. Л.: Гидрометеоиздат, 1987. 254 с.
[3] Maltsev V.P. Scanning flow cytometry for individual particle analysis // Review of Scientific Instruments. 2000. Vol. 74. P. 243-255.
[4] Апонасенко А.Д., Франк Н.А., Сидько Ф.Я. Дифференциальный спектрофотометр для гидрооптических исследований // Океанология. 1976. Т. 16, № 5. С. 924-928.
[5] Hunt A.J., Huffman D.R. A new polarization-modulated light scattering instrument // Review of Scientific Instruments. 1973. Vol. 44. P. 1733-1762.
[6] Варшалович Д.А., Москалев А.Н., Херсонский В.К. Квантовая теория углового момента. Л.: Наука, 1975. 439 c.
[7] Kuscer I., Ribaric M. Matrix formalism in the theory of diffusion of light // Opt. Acta. 1959. Vol. 6. P. 42-51.
[8] Hovenier J.W., van der Mee C.V.M. Fundamental relationships relevant to the transfer of polarized light in a scattering atmosphere // Astron. Astrophys. 1983. Vol. 128. P. 1-16.
[9] Mishchenko M.I., Travis L.D., Lacis A.A. Scattering, Absorption, and Emission of Light by Small Particles. Cambridge: Cambridge Univ. Press, 2002. P. 445.
[10] Гельфанд И.М., Минлос Р.А., Шапиро З.Я. Представления группы вращений и группы Лоренца и их применения. М.: ГИТТЛ, 1958. 368 с.
[11] Paramonov L.E. T-matrix approach and the angular momentum theory in light scattering problems by ensembles of arbitrarily shaped particles //J. Opt. Soc. Am. A. 1995. Vol. 13. P. 2698-2707.
[12] Paramonov L.E. Light scattering by randomly oriented particles into solid angles //J. Opt. Soc. Am. A. 1994. Vol. 11, N 4. P. 1360-1369.
[13] Waterman P.C. Symmetry, unitarity and geometry in electromagnetic scattering // Phys. Rev. D. 1971. Vol. 3, N 4. P. 825-839.
[14] Tsang L., Kong J.A., Shin R.T. Radiative transfer theory for active remote sensing of layer of nonspherical particles // Radio Sci. 1984. Vol. 19, N 2. P. 629-642.
[15] Абдулкин В.В., Парамонов Л.Е. Решения векторного уравнения Гельмгольца, инвариантные относительно группы вращений // Вопр. мат. анализа: Сб. науч. тр. / Красноярский гос. техн. ун-т, 2004. Вып. 7. С. 3-15.
[16] Шмидт В.А., Парамонов Л.Е. Фурье-разложение элементов матрицы рассеяния хаотически ориентированных несферических частиц // Вопр. мат. анализа: Сб. науч. тр. / Красноярский гос. техн. ун-т, 2004. Вып. 7. С. 154-164.
[17] Борен К., Хафмен Д. Поглощение и рассеяние света малыми частицами. М.: Мир, 1986. 660 c.
[18] Mishchenko M.I. Light scattering by randomly oriented axially symmetric particles //J. Opt. Soc. Amer. A. 1990. Vol. 8. P. 871-882.
[19] Гуревич М.М. Фотометрия (Теория, методы и приборы). Л.: Энергоатомиздат, 1983. 268 с.
[20] Парамонов Л.Е., Осколкова Г.В. Рассеяние света в телесных углах при сфокусированном или расходящемся падающем излучении // Оптика и спектроскопия. 1993. Т. 74, вып. 2. С. 176-182.
Поступила в редакцию 2 ноября 2005 г.
