Глава 2. Метод дискретных диполей

2.1.Обзор МДД*

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

2.1.1. Введение

МДД это метод моделирования рассеяния и поглощения электромагнитных волн частицами произвольной формы и внутренней структуры. Изначально МДД был предложен Пурселлом (Purcell) и Пеннипэкером (Pennypacker) (ПП) [44], которые исходили из физической картины набора точечных диполей. В дальнейшем он развивался Дрейном (Draine) и сотрудниками [45-48], которые пропагандировали метод путём развития общедоступной компьютерной программы DDSCAT (см. параграф 2.5.2.2). Позднее было показано, как вывести МДД из интегрального уравнения для электрического поля путём его дискретизации разделением всего объёма рассеивателя на малые кубические элементы. Этот вывод был впервые проведён, вероятно, Годеке (Goedecke) и О’Браеном (O’Brien) [49] и в дальнейшем развивался другими [50-53]. Важно отметить, что оба вывода МДД приводят к одним и тем же уравнениям, различаясь лишь в том, что формулировка на основе интегральных уравнения даёт больше математического понимания используемых приближений, тем самым показывая возможные пути улучшения метода, а модель точечных диполей более понятна физически.

Некоторые исследователи называют МДД методом (или приближением) связанных диполей (coupled dipole method) [54,55]. Другие методы, такие как формулировка на основе объёмного интегрального уравнения (volume integral equation formulation) [56] и дискретизированная функция Грина (ДФГ, discretized Green’s function) [49], разрабатывались независимо от ПП, однако позднее была показана их эквивалентность МДД [27,50]. В дальнейшем, термин МДД будет использоваться для всех этих методов, так как они попадают под общую формулировку. Тем не менее

* Данный обзор опубликован в Yurkin M.A., Hoekstra A.G. The discrete dipole approximation: an overview and recent developments. // J. Quant. Spectrosc. Radiat. Transf. – 2007. – V.106 . – P.558-589.

19

тяжело однозначно отделить МДД от методов основанных на объёмном интегральном уравнении, таких как широкий диапазон ММ с различными базисными и тестовыми функциями [57-60]. Мы считаем, что фундаментально МДД отличается тем, что неизвестные величины, определение которых является первым и главным этапом в процессе решения, являются внутренним электрическим полем или его производной величиной, например, поляризацией, что поддаётся прямой физической интерпретации. Иными словами, любая разновидность МДД может быть интерпретирована как замена рассеивателя множеством взаимодействующих диполей, что обсуждается в подразделе 2.1.2. В качестве примера метода, который мы не относим к МДД, можно привести ММ с иерархическими базисными функциями Лежандра высокого порядка [58].

МДД известен в области светорассеяния и он описан в нескольких обзорах. Обширный обзор Дрейна и Флато (Flatau) [47] описывает практически всё развитие МДД до 1994 г. Более современный обзор Дрейна [48] в основном направлен на применения и численные аспекты. Теория МДД, совместно с другими методами, рассмотрена в обзорах Рида (Wriedt) [61], Чиапетты (Chiappetta) и Торезани (Torresani)

[62]и Канерта [27], а также в монографиях Мищенко и др. [26] и Цанга (Tsang) и др.

[63].Джонс (Jones) [64] рассматривал МДД в контексте различных методов для характеризации частиц. Тем не менее многие улучшения МДД предложенные после 1994 г. не упомянуты ни в одной из этих работ. Те же, что упомянуты, обычно считаются боковыми отступлениями и не рассматриваются в рамках общей формулировки. Более того, насколько нам известно, численная сторона МДД нигде не описана подробно – в каждой работе рассматриваются лишь несколько конкретных аспектов. В данном обзоре мы постараемся заполнить эти пробелы.

Общая формулировка описывается в подразделе 2.1.2 для облегчения дальнейшего обсуждения различных вариантов МДД. Эта формулировка основана на интегральном уравнении, что позволяет единообразно описать всё развитие МДД, однако связь с физически понятной моделью точечных диполей также чётко обозначена. Там же обсуждаются источники ошибок в МДД.

Вподразделе 2.1.3 описаны физические основы МДД и сравнены результаты различных вариаций метода. В параграфе 2.1.3.1 с теоретической точки зрения описаны различные выражения для поляризуемости, члена, описывающего взаимодействие, и сечения поглощения Cabs. Сравнение результатов моделирования, используя различные выражения, приведено в параграфе 2.1.3.2. В частном случае кластеров шаров возможны некоторые улучшения и упрощения, которые описаны в параграфе 2.1.3.3, а существенные модификации, выходящие за рамки общей формулировки, описанной в

20

подразделе 2.1.2, или предназначенные для узких целей, обсуждаются в параграфе

2.1.3.4.

В подразделе 2.1.4 рассматриваются различные численные аспекты МДД, связанные, в основном, с решением очень больших систем линейных уравнений (параграф 2.1.4.1). Параграф 2.1.4.2 посвящён простому итерационному методу для решения линейной системы МДД, обладающему ясным физическим смыслом. Специальная структура матрицы взаимодействия при использовании прямоугольной решётки и применение этой структуры для ускорения вычислений обсуждается в параграфах 2.1.4.3 и 2.1.4.4 соответственно. Общие методы для ускорения вычислений, не требующие прямоугольной решётки, рассматриваются в параграфе 2.1.4.5, а параграф 2.1.4.6 посвящён специальным техникам увеличения эффективности повторяющихся вычислений (например, при усреднении по ориентации).

Численные сравнения МДД с другими методами рассмотрены в подразделе 2.1.5; там же обсуждаются его сильные и слабые стороны. подраздел 2.1.6 завершает обзор и представляет возможное дальнейшее развитие МДД.

2.1.2. Общая формулировка

Везде в данной диссертации предполагается что временная зависимость всех полей описывается множителем exp(−iωt), а рассеиватель обладает диэлектрическими, но не магнитными свойствами (единичная магнитная проницаемость). Для упрощения выкладок предполагается, что диэлектрическая проницаемость изотропна, но рассмотрение можно обобщить на случай произвольного тензора диэлектрической проницаемости.*

Электрическое поле внутри диэлектрического рассеивателя описывается следующим интегральным уравнением [27,50]:

χ(r) = (ε(r) − 1)/4π – восприимчивость среды в точке r (ε(r) – относительная проницаемость). V обозначает объём частицы, т.е. объём, состоящий из всех точек, где восприимчивость ненулевая, а V0 – меньший объём такой, что V0 V, r V0\∂V0 (∂

* В большинстве формул скалярные величины можно напрямую заменить тензорным, но есть исключения. Расширения МДД для оптически анизотропных рассеивателей обсуждается в параграфе

2.1.3.4.

21

обозначает границу области). G(r,r′) обозначает тензорную функцию Грина в свободном пространстве, определяемую как

M обозначает следующий интеграл, связанный с конечностью исключаемого объёма V0

где nˆ′ обозначает внешнюю нормаль к поверхности ∂V0 в точке r'. Тензор L всегда вещественный и симметричный, а его след равен 4π [65]. Важно отметить, что L не зависит от размера V0, а только от его формы (и положения точки r внутри), в то время как M зависит от размера элемента объёма и стремится к нулю при его уменьшении

[50] (если и χ(r), и E(r) непрерывны внутри V0).

При выводе уравнения (2) особенность функции Грина выделяется в явном виде, поэтому оно предпочтительнее чем часто используемое выражение [27,50]:

V

Более того, Янгджан (Yanghjian) отметил [65], что существует несколько способов рассмотрения особенности в уравнении (8), приводящих к разных результатам. Он также доказал, что вывод уравнения (8) некорректен в окрестности особенности функции G(r,r′) . Поэтому, оно может считаться верным, только если особенность раскрывается так же, как это делал Лахтакиа (Lakhtakia) [50], сводя его к строгому уравнению (2).

22

N

Уравнение (2) дискретизируется следующим образом [27]: пусть V = UVi ,

i=1

Vi IVj = 0/ для i ≠ j, где N обозначает количество объёмных элементов.* Хотя в данной

формулировке объёмные элементы могут быть любыми, в большинстве случаев используются одинаковые стандартные элементы (ячейки). Обычно набор этих ячеек не может точно описать форму рассеивателя, и дискретизация является приближённой.

Предполагая r Vi и выбирая V0 = Vi, переписываем уравнение (2) в виде

Система уравнений (9) (для всех i) является точной. Далее выбираем по одной точке ri в

каждом Vi (его центр) и устанавливаем r = ri. Во многих случаях можно сделать следующие предположения:

Vj

(11)

которые состоят в том, что интегралы в формуле (9) линейно зависят от значений χ и E в точке ri. Тогда уравнение (9) переписывается в виде

Верхний индекс (0) обозначает приближённое значение тензора (при некоторых предположениях). Практически все разновидности МДД, например [50], используют ещё одно приближение:

* В данной формулировке МДД объёмный элемент также называется диполем.

23

Это предположение неявно используется всеми видами МДД, которые начинают с замены рассеивателя множеством точечных диполей. Важно отметить, что уравнение (12) и выводы из него требуют более слабых предположений (10) и (11) чем обычно используемые предположения (13) и, тем более, (16). Например, формулировка Пелтониеми (Peltoniemi) [66], описанная в параграфе 2.1.3.1 основана непосредственно на уравнении (12). Мы постулируем уравнение (12) как отличительную черту МДД, т.е. метод относится к МДД тогда и только тогда, когда его главное уравнение эквивалентно уравнению (12) с некоторыми Vi, χi, Mi , Li и Gij .

Канерт [27] разделил МДД и ММ на основании того, что ММ непосредственно решает уравнение (12), находя неизвестные Ei, в то время как МДД находит не суммарное, а возбуждающее электрическое поле

где Eselfi обозначает самонаведённое поле объёмного элемента. Тогда уравнение (12)

эквивалентно

Существует, однако, альтернативная формулировка МДД [47], которая находит неизвестную поляризацию диполей Pi:

(21)

(22)

Важно отметить, что Pi, определяемая формулой (21), лишь приближение для полной поляризации элемента Vi, которое точнó только при условии (13), в то время как формулировка в целом этого условия не требует. Можно рассматривать уравнение (22) как промежуточное звено между МДД и ММ в классификации Канерта [27], что показывает полную эквивалентность всех этих формулировок. Уравнение (22) предпочтительней уравнений (12) и (19) для численного решения, ввиду специальной структуры матрицы Gij , что подробно описано в подразделе 2.1.4.

Лахтакиа [50] определил «сильную» и «слабую» формулировки МДД как те, что соответственно учитывают или пренебрегают Mi . Эти два варианта сближаются при

24

Совместно с формулой (20) это приводит к известной формулы для поляризуемости Клаузиуса-Моссотти (KM) для слабой формулировки МДД:

где εi = ε (ri), а d обозначает размер кубической ячейки. Эта же поляризуемость использовалась Пурселлом и Пеннипэкером [44].

После вычисления внутреннего поля, можно посчитать рассеянное поле и интегральные сечения. Рассеянное поле получается в пределе r → ∞ из интеграла в формуле (2) (см., например, [49]):

где n = r/r это единичный вектор в направлении рассеяния, а F – амплитуда рассеяния:

i Vi

Все другие дифференциальные характеристики рассеяния, такие как амплитудная матрица рассеяния и матрица Мюллера, а также параметр асимметрии <cosθ >, можно получить из F(n), посчитанного для двух поляризаций падающей волны [14]. Можно также вычислить силу, с которой излучение действует на рассеиватель целиком или на его части [36,67,68]. Рассмотрим падающее излучение в виде плоской волны*

Einc (r) = e0 exp(ik r) , (27)

где k = ka, a – направление распространения, а |e0| = 1. Тогда сечение рассеяние определяется [14] как

Сечения поглощения и экстинкции (Cabs, Cext) вычисляются [49,56] непосредственно из внутреннего поля:

iVi

*МДД подходит для любого падающего поля, например для гауссового пучка [69]; но здесь это не обсуждается.

25

Однако, как было отмечено Дрейном [45], при вычисления Csca формула (31) может приводить к бóльшим погрешностям чем (28), особенно при Cabs >>Csca.

Самый простой способ представить выражения (26) и (29) в терминах внутреннего

Оставляя только основной член в разложении экспоненты около ri, можно упростить формулу (32) до

(34)

что, совместно с формулой (30), приводит к

Cext(0) = 4πk∑Im(Pi Einci ). (35)

Формулы (34) и (35) идентичны тем, что использовались Пурселлом и Пеннипэкером [44], а затем и Дрейном [45], выражения же для Cabs слегка отличаются, что обсуждается в параграфе 2.1.3.1. К сожалению, многие исследователи явно не указывают, как именно они вычисляют характеристики рассеяния из внутреннего поля или поляризации диполей. Те, что указывают, обычно используют процедуру Дрейна

[формулы (28), (34), (35) и (37)].

Погрешности МДД могут быть разделены на ошибки дискретизации, связанные с конечностью размера ячейки d, и ошибки формы, связанные с неточностью описания формы частицы набором стандартных ячеек, например, кубических. Ошибки дискретизации возникают в результате предположения постоянства E внутри каждой ячейки и приближённого вычисления Mi и Gij . Можно считать, что ошибки формы следуют из предположения постоянства χ и E внутри граничных ячеек, что неверно, так как эти ячейки пересекаются поверхностью частицы. С другой стороны, можно рассматривать ошибку формы как разницу между точными результатами для частицы исходной формы и для частицы, составленной из набора стандартных ячеек. Оба типа ошибок стремятся к нулю, когда N → ∞ при постоянных геометрии рассеивателя и

26