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

УДК 681.325.5:681.513.6

А.М. Серебренников

Пермский национальный исследовательский политехнический университет

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ МУЛЬТИПОЛЬНЫХ ВЗАИМОДЕЙСТВИЙ В ЦЕПОЧКАХ ЧАСТИЦ С ПОМОЩЬЮ ОДНОГО ЧИСЛЕННО-АНАЛИТИЧЕСКОГО МЕТОДА

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

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

в которой будет распространяться информация, и при этом получать очень компактные устройства. Конструктивно такие среды могут быть одно-, двух- и трехмерными. Ожидается, что в этой области появятся новые идеи по части гибридизации электронных устройств и искусственных материальных сред для передачи и обработки информации. Разумеется, потребности техники диктуют необходимость создания новых средств разработки (компьютерных электромагнитных симуляторов), в основе которых лежали бы эффективные математические модели, методы и алгоритмы. Поскольку подходящим элементарным строительным блоком для создания плазмонных устройств рассматриваются сферические частицы, представляется разумным развивать методы, базирующиеся на сферическом анализе. Мы не говорим здесь о применении теории Ми в чистом виде [1], поскольку в любых реалистичных технологиях существуют производственные допуски и любые устройства имеют неконтролируемые геометрические отклонения в пределах этих допусков. Влияние таких отклонений на функциональность устройств необходимо принимать во внимание. В результате требуется моделировать ситуации с частицами несферической (а лучше произвольной) формы. В статьях [2-4] мы рассматривали различные подходы к решению задач рассеяния (излучения) с диэлектрическими частицами. В [3] нам удалось создать достаточно универсальные алгоритмы для работы с частицами сложной геометрии. В этой статье мы обобщим одну из наиболее удачных систем уравнений, полученную для одночастичной задачи в [3], на случай системы многих тел. Покажем, что учет трансляционной симметрии приводит к особому способу формирования разрешающей системы уравнений и позволяет получить существенный вычислительный выигрыш. Также мы решим некоторые тестовые задачи и проведем внутренние проверки в обеспечение верификации наших методов и алгоритмов. Кроме того, в работе проводится обсуждение разработанного нами метода итерационного решения переопределенных систем линейных алгебраических уравнений (СЛАУ). Приводится обоснование его сходимости. Наконец, мы сфокусируемся на применении разработанных моделей и методов к исследованию цепочечных систем, представляющих интерес в качестве оптических волноводов [5-7]. Их исследованию посвящена секция численных результатов.

Теория, постановка задачи и метод решения. Поскольку в работе мы остаемся в рамках линейной электродинамики изотропных сред, то постановка задачи сохраняется в основном той же, что и в предыдущих статьях [2-4]. В целях сокращения текста она не приводится. Добавляются два новых существенных пункта. Во-первых, новую специфику в формулировку задачи вносит физика явления, моделирование которого производится в данной статье. А именно колебания газа коллективизированных электронов в металле под действием периодического поля световой волны являются чисто квантовым явлением и объясняются на основе квантово-механических теорий в физике конденсированного состояния вещества. Однако за более чем столетнюю историю изучения этого явления было замечено, что электронные состояния в металлах в области оптических частот хорошо описываются моделью газа свободных электронов (моделью Друде).

Наиболее просто эта модель может быть получена на основе второго закона Ньютона [8], но также является частным случаем и из других физических теорий. В частности, она следует из феноменологической теории Ландау [9], а также из микроскопических теорий [10], [11]. В соответствии с моделью Друде диэлектрическая функция плазмонного металла имеет вид:

ег (ю) = 1 - юр [ю (ю + Ще1 )]-1,

где юр1 - плазменная частота, юге1 - параметр релаксации. Эти параметры следует понимать как принципиальные материальные константы теории, значения которых добываются экспериментально. К настоящему моменту они известны для основных металлов, обладающих плазмоподобным поведением (А§, Аи, Си). Во-вторых, различие в формулировке обусловлено тем, что сейчас мы решаем задачу многих тел. По этой причине мы воспользуемся формализмом линейных волновых теорий.

Пусть наша система состоит из N частиц (тел). С каждой частицей мы свяжем ее собственную локальную «волновую функцию» (ВФ) в виде мультипольного ряда или (что эквивалентно) ряда Аткинсона-'^1сох [2], [12]. Такой подход к решению вполне очевиден, он используется и в других численных методах электродинамики [13]. Параллельно на нем базируются и некоторые методы теории конденсированного состояния вещества. Отметим среди них метод ЛКАО [14],

применяемый в квантово-химических численных моделях для расчета электронных состояний молекул. (В физике кристаллов этот метод известен также как метод сильной связи [15].) В нем одночастичная ВФ электрона сложной многоатомной системы строится в виде линейной комбинации примитивных ВФ - атомных орбиталей.

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

первичное поле Е0,Н0. Далее это поле взаимодействует с системой

тел, в результате чего во вмещающей среде индуцируется рассеянное

Т^БО ^Ш^N т-’вХ ТТ БО '^~УN ТТ еХ

поле Е = Ед , Н = ^=1 Нч , где под знаком суммы стоят одночастичные поля Ее*, Не* (одночастичные ВФ). Внутри частиц воз-

ттш

никают в свою очередь внутренние поля Ер , Н , которые принципиально «одночастичны». С учетом этих обозначений граничные условия на поверхностях частиц Бр можно записать как

п х (Е0 + Е30) = п х Ер , п х (Н0 + Н30) = п х Нр , (1)

' 0 ' п рС ' 0 ' ~ рс

Бр Б р Бр Бр

где п обозначает единичную нормаль к Бр, направленную из частицы

во внешнюю среду. Для построения разрешающей системы уравнений мы так же, как и в [3], используем интегралы Стрэттона-Чу [12], описывающие внутренние поля в частицах. Возьмем произвольную частицу (с индексомр) и запишем для нее интеграл внутреннего поля в ДСК:

I н р {г')т йу=

<2р

(2)

= №Р х (п х и;) йБ + -— jgppp (п х Ер) йБ,

бр кр2р бр

где Я ° | г - г' |, gp = ехр(ikpR)/4л Я - функция Грина однородного пространства, 3(Я) - дельта-функция Дирака,

Рр = - Я-2 (1 - ikpR - к2р Я2) и - Я-2 (-3 + ШрЯ + к2р Я2) еяея ,

где и - единичный тензор, Я ° | г - г' |, еЯ ° (г - Г)/ Я, еЯеЯ - диад-ное произведение, г - точка интегрирования, Г - точка наблюдения поля, кр - волновое число; 2р - волновой импеданс. Подставим (1)

в (2) и поместим точку наблюдения г' за пределы объема Qp, занятого частицей, во внешнюю область. В этом случае благодаря свойствам дельта-функции объемный интеграл в левой части обратится в нуль.

При представлении полей Е“, Ие* в виде рядов Аткинсона-Wilcox уравнение (2) приобретет вид интегрального:

I% ххй")+ 1gpPp хЕе)=

Бр р р Б

= х (п хи0) |gpрр(п хЕ0) <Ж.

Бр р р Бр

р р Бр (3)

ррр(п хЕ0,)

При представлении Ее*, Ие* мультипольными разложениями

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

Е;к,] = к2[у9]+ikzKl[uq], Н;х[и,V] = К2"]-ikz-lК^], (4)

удовлетворяющих однородным уравнениям Гельмгольца, называются векторным мультипольным (МП) разложением ЭМ поля, если скалярные поля и и vq выражены посредством рядов вида

\"; ^ = I Ь"" } ^ (вq , Ф, ), (5)

[ ^ V )пт пт

где к!1 - сферическая функция Ханкеля порядка п,

¥пт , Фq ) = ^ ) - сферическая фуHкЦия, Д|т| (сов eq ) -

нормированный присоединенный полином Лежандра, К1[ ] °Ух(г )

и К2[ ] °УхК1 [ ] - дифференциальные операторы в локальной сферической системе координат. Суммирование в (5) выполняется по (п,т), при этом п = 0,...,¥ и т = -п,...,п . Если выразить более точно, ряды (4)

являются разложением внешних полей Еечх, Ие*. Внутренние поля Ер,

Ир выражаются аналогичной зависимостью с той лишь разницей, что

h:\kr,) меняется на ]п(кргр), т.е. сферическую функцию Бесселя,

а материальные константы без индексов, относящиеся к внешней среде, - на соответствующие константы среды, заполняющей частицы. (Индексы р и q в обозначениях расстояний и углов показывают их принадлежность к локальным системам координат, связанным с частицами, изображенными на рис. 1).

Поскольку положение точки г' в уравнении (3) может быть выбрано произвольным (с сохранением требования: г'£ Qp

и г' £ Бр ), то, размещая г' в различных точках на некоторой вспомогательной (виртуальной) поверхности, описанной вокруг р-й частицы, мы наберем систему линейно-независимых алгебраических уравнений относительно МП коэффициентов. Линейная независимость относительно коэффициентов " и vqnm следует из того, что координаты точки г входят в уравнение (3) нелинейно. Все остальное, что можно сказать по поводу уравнения (3), вполне аналогично тому, что обсуждалось в [3]. Таким образом, уравнения (3), собранные вме-

сте, образуют переопределенную алгебраическую систему вида МЖ = V или более подробно:

(6)

"Мц . •• К • •• Л~ті с -*■ л VI

Мрі • - ^ • « % > = — « к

_Мт • •• Мтч • • ■ ^~тт Vт

где все элементарные матрицы М имеют большее число строк, чем

столбцов. В основе численного метода решения системы (6) лежит процедура минимизации квадратичной невязки:

8(Ж )

V - МЖ

—*гр —*■ ~*гр г~ —*■ —■► гг, гг, —*■

VTV - 2 VіМ Ж + ЖТМТМ Ж, (7)

где Т обозначает транспонирование вместе с комплексным сопряжением, ||*|| л - норму комплексного арифметического векторного евклидова пространства, с1 - размерность вектора V . Условие минимума ( дЪ/дЖ = 0) ведет к нормальной системе:

ММТМ Ж = ММTV, (8)

которая может быть решена любым численным методом. Решение нормальной системы доставляет минимум невязке (7):

8ШП = VTV - VіМ Ж = VTV - ЖТМТМЖ.

(9)

За счет увеличения порядка конечной суммы МП ряда величина Ътт может быть минимизирована до сколь угодно малого (заданного наперед) уровня. Доказательство этого утверждения приводится в [3]. Для этого требуется некоторое преобразование базисных функций и коэффициентов МП ряда. Техника этого преобразования универсальна и не изменяется при переходе от одночастичной задачи к задаче многих тел. При разработке алгоритма решения системы (6) заслуживает внимания лишь одна техническая деталь. Дело в том, что МП ряды (4), (5) записаны в локальных сферических системах координат, ассоциированных с частицами (см. рис. 1). Таким образом, при связывании их посредством уравнения (3) требуется процедура преобразования компонент векторных мультиполей в единую (декартову) систему координат.

2

а

С

Технически это делается путем умножения векторов одночастичных полей Е“, Не на некоторые матрицы преобразования. Технику таких вычислений мы здесь опускаем, так как она хорошо известна в векторном (тензорном) анализе и может быть легко найдена в соответствующей литературе (см., например, [16]).

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

удовлетворяющих условию Я ^ = Ярч, означает появление трансляционной симметрии. Это, в свою очередь, означает тождественность матриц М^ и Мп . Эта особенность является исключительно важной

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

большое число совпадений вида = Ярд = Яи = • • • = Я^. Таким образом, в алгоритме решения требуется вычислять лишь неповторяющиеся блоки и располагать их в глобальной матрице СЛАУ (6) на соответствующих позициях. Эффективность учета трансляционной симметрии мы продемонстрируем на примере прямой регулярной цепочки, состоящей из N частиц. При разупорядочении цепочки количество различных элементарных матриц, которые необходимо сформировать, очевидно, будет равно квадрату N. Непосредственным подсчетом легко установить, что в регулярном случае число различных матриц в системе (6) составляет всего 2N -1 штук. Благодаря этому мы можем говорить о линейном алгоритмическом масштабировании (скейлинге) вместо квадратичного. Конкретная практика вычислений подтвердила высокую эффективность учета трансляционной симметрии. Однако в случае больших систем уравнений мы получаем ощутимый выигрыш только на этапе формирования СЛАУ. Строго говоря, нам еще требуется эффективный метод решения алгебраической системы. Его разработке посвящен следующий раздел.

Итерационный метод решения переопределенных СЛАУ. С вычислительной точки зрения у переопределенных систем есть один существенный недостаток. Для решения переопределенной СЛАУ прямым методом ее необходимо предварительно привести к нормальному виду (8). Это достаточно затратная по времени проце-

дура, требующая порядка ~ ш2N операций умножения только для вычисления матрицы МТМ (где N - число строк, а ш - число столбцов исходной матрицы М). Разумеется, после этого прямой метод решения нормальной системы (например, метод Гаусса) добавляет

еще ~ ш3 операций. Процедура умножения Мт на М происходит достаточно быстро, если ш << N, но для матриц, близких к квадратным, алгоритм оказывается слишком медленным. Было бы желательно найти какой-нибудь итерационный подход, в котором было бы не очень много операций умножения в пересчете на одну итерацию.

Предлагается следующий. Возьмем квадратичную невязку 8(Ж) и дадим приращение ее аргументу:

2

~»а

Ъ(Ж + ёЖ) = (V - МЖ) -МёЖ = Ъ(Ж) - 2 ёЖТМТ (V - М Ж) + ёЖТМТМёЖ.

(Эту функцию можно понимать как точное разложение Ъ(Ж) в степенной ряд.) Отрицательность конечной разности Ъ(Ж + ёЖ) - Ъ(Ж) при заданном ёЖ гарантирует движение в сторону минимума Ъ(Ж) . Таким образом, реализуется итерационная схема вида Жг+1 = Жг + ёЖг,

где г - итерационный шаг алгоритма. Неположительность конечной разности требует на каждом шаге выполнения условия:

- 2 ёЖТМТ (V - МЖ) + ёЖТМТМёЖ < 0. (10)

Пусть ёЖ = 2М (V - МЖ)/С , где С > 0 - некоторая константа, подлежащая определению. (Очевидно, что при такой записи вектор ёЖ представляет собой нормированный на С вектор антиградиента.) С учетом принятых обозначений неравенство (10) принимает вид: - С || ёЖ ||2 + || МёЖ ||2 < 0, откуда мы получаем условие:

С > || МёЖ ||2 || ёЖ ||-2. (Размерность пространства в обозначении нормы опускаем.) Применяя неравенство Коши-Буняковского || МёЖ ||2 < || М ||21| ёЖ ||2 и решая совместно два последних неравенства, получаем: С > || М ||2 . Данное неравенство следует понимать как условие сходимости метода. Итерационная процедура нахождения

минимума невязки 8(Ж) на первом шаге работы алгоритма может быть организована как

Ж = Ж + ёЖ = Ж + 2МТ (V- МЖ) а-1 ||М~||-2,

где а > 1, а Ж - начальное приближение. На последующих шагах модифицируем ее как

—*• ^ _ —*• ^ —*■ 1 —*• _ —*• _

Ж+1 = Ж + ёЖг = Ж + 2М (V - М^) а-1 ЦёЖ^ХЦМёЖ,-1||-2.

Критерием остановки алгоритма является минимизация 8(Ж)

(а также ёЖ ) ниже некоторого наперед заданного уровня. Практика использования данного алгоритма показала необходимость формулировки метода генерации подходящего начального приближения. В противном случае количество итераций оказывается достаточно большим. Ясно, что для произвольного случая создать такой метод не представляется возможным. Однако для волноводных задач с цепочками частиц, облучаемых компактными источниками, оказалось возможным создать достаточно простой и эффективный подход. Предлагается следующий эвристический итерационный метод. (Принципиальный момент состоит в том, что его сходимость не может быть доказана в общем случае.) Идея метода базируется на физических соображениях. Для определенности рассмотрим задачу о возбуждении продольного плазмон-поляритона вдоль однорядной цепочки частиц компактным источником, находящимся в окрестности одной из частиц подобно тому, как это показано на рис. 2.

Источник

2 I

Рис. 2. Продольная плазмонная мода

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

пени взаимодействием с соседней частицей. Матрица первой частицы М11 занимает левый верхний угол глобальной матрицы М . Алгоритм

строится так. На первом шаге решаем систему: М11 #1 = Ух. Вторая частица взаимодействует с полем источника и с первой частицей, поле которой теперь уже найдено, третья частица взаимодействует с частицами 2 и 1 и т.д. В результате строится следующий алгоритм последовательного нахождения одночастичных векторов-решений:

М~22 ^2 = V2 — М21 ^1,

М33 Щ = У3 — А~31 ^1 — М~32 ^2 ,

М„™, = V,— 'Е—=1М .

После того как все элементарные векторы определены,

строится следующий итерационный процесс:

ДЛ- = V,— ЕтФІМу™і.

Критерием его остановки вновь является минимизация 8(Ж) (а также ёЖ) вплоть до минимально достижимого уровня. Далее включается алгоритм градиентного поиска минимума невязки 8(Ж), описанный выше.

ЭПР, дБ

18 -I--------т---------,---------,---------,--------,----------------

О 30 60 90 120 150 180 0

Рис. 3. Диаграмма ЭПР

Численные результаты. Эту секцию мы начнем с решения одной тестовой задачи в обеспечение верификации нашего метода и алгоритмов. Рассмотрим диэлектрический куб с параметрами Ь = 1,0 м

(длина ребра), ег2 = 3,0, тг2 = 1,0, = 0 см/м, облучаемый плоской

волной частотой 299,792 МГц. Волна предполагается падающей нормально из вакуума на одну из граней куба. На рис. 3 показана диаграмма эффективной поверхности рассеяния (ЭПР), полученная несколькими альтернативными численными методами.

Кривые 1 и 2 получены с помощью систем уравнения (9) и (27) из статьи [3], кривая 3 - с помощью уравнения (31) из [2], кривая 4 -с помощью стандартного метода Т-матриц [17]. Далее мы соберем этот же куб из восьми кубов меньшего размера (Ь = 0,5 м) и решим для такой сборки задачу многих тел. На рис. 3 кривая 5 получена с помощью уравнения (31) из работы [2], обобщенного на случай системы многих тел, а кривая 6 - с помощью системы (6) данной статьи. Видно, что между всеми кривыми нет визуальных различий. Сейчас рассмотрим задачу о распространении продольной волны плазменных колебаний вдоль однорядной цепочки, состоящей из десяти сферических малых частиц (Б = 50 нм), показанной на рис. 2. Расстояния между центрами частиц Ь = 55 нм. Такая волна может быть создана источником, генерирующим электрическое поле, направленное по оси Х. В качестве источника возьмем отрезок тока, расположенный рядом с поверхностью крайней левой частицы. На рис. 2 показано поверхностное распределение напряженности электрического поля 1§ | Ё*с |.

На рис. 4 показано распределение напряженности в центрах частиц и зазорах между ними. С помощью аналитического интегрирования системы (31) из [2] можно показать, что для малых частиц коэффициенты внутреннего МП разложения есть апт ~ г—”+1, Ьпт ~ г—”+1 (г - радиус частицы), и при этом они нелинейно зависят от ю, &р1, . Таким образом, резонансному усилению у нас подвергается почти полная группа мультиполей (п > 2 ) и при этом неоднородно.

Эта особенность объясняет сложный характер резонансного поведения наночастиц (плазмонный резонанс), известный по многочис-

ленным исследованиям (см., например, [13]). Явление сильной концентрации мы наблюдаем и в нашем примере. Быстрый спад напряженности, как видно из рисунков, является следствием омических потерь, так как поле «переконцентрировано» в частицах вследствие плазмонного резонанса и межчастичных взаимодействий, а релаксационный параметр юге1 не является достаточно малым. Параметры модели Друде соответствуют серебру и обобщены по результатам анализа множества литературных источников, в том числе тех, что указаны в списке литературы в конце статьи. В расчетах использовались следующие значе-

ния: юр1 = 2,5 -1015 Гц, ю, = 2,0 1013 Гц, ю = 4,5 1014 Гц.

,13

14

р1

К

н

те

ч

о

с

н

о

о

к

к

(и

*

те

ср

с

Й

Я

Номер частицы

Рис. 4. Электрическое поле в волноводе: 1 - напряженность электрического поля в центрах частиц для волноводной моды, изображенной на рис. 2;

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

В таблице приводятся результаты исследования сходимости метода. В первой колонке помещены значения старшей степени волновой гармоники в МП разложении, во второй - нормированной невязки (9), в третьей - полной «дальнезонной» излученной мощности, в четвертой - напряженности электрического поля в точке А (см. рис. 2), т.е. в зазоре между частицами.

Результаты исследования сходимости метода

Nmultip 8 . • || V |Г2 min II ' II Pscat , отн. ед. Ea , отн. ед.

7 5.3Е-8 1795.3 2443.6

8 8.0Е-9 1787.2 2437.5

9 9.0Е-9 1789.8 2465.3

10 9.0Е-9 1788.0 2461.2

Библиографический список

1. Борен К., Хафмен Д. Поглощение и рассеяние света малыми частицами. - М.: Мир, 1986.

2. Serebrennikov A.M. An analysis of scattering caused by dielectric bodies using semi-analytic methods // IEEE Trans. Antennas Propagat. -2008. - Vol. 56. - P. 3201-3209.

3. Serebrennikov A.M. A novel semi-analytic method for the analysis of scattering by dielectric objects immersed in uniform media // Comp. Phys. Comm. - 2010. - Vol. 181. - P. 1087-1095.

4. Серебренников А.М. К обоснованию сходимости одного мультипольного метода для решения векторных задач рассеяния. Системы мониторинга и управления: сб. науч. тр. - Пермь: Изд-во Перм. гос. тех. ун-та, 2010. - С. 203-207.

5. Park S.Y., Stroud D. Surface-plasmon dispersion relations in chains of metallic nanoparticles: An exact quasistatic calculation // Phys. Rev. B. - 2004. - Vol. 69.

6. Weber W.H., Ford G.W. Propagation of optical excitations by dipolar interactions in metal nanoparticle chains // Phys. Rev. B. - 2004. - Vol. 70.

7. Govyadinov A.A., Markel V.A. From slow to superluminal propagation: Dispersive properties of surface plasmon polaritons in linear chains of metallic nanospheroids // Phys. Rev. B. - 2008. - Vol. 78.

8. Слэтер Дж. Диэлектрики, полупроводники, металлы. - М.: Мир, 1969.

9. Бредов М.М., Румянцев В.В., Топтыгин И.Н. Классическая электродинамика. - СПб.: Лань, 2003.

10. Шриффер Дж. Теория сверхпроводимости. - М.: Наука, 1970.

11. Маделунг О. Теория твердого тела. - М.: Наука, 1980.

12. Колтон Д., Кресс Р. Методы интегральных уравнений в теории рассеяния. - М.: Мир, 1987.

13. Stout B., Auger J.C., Devilez A. Recursive T matrix algorithm for resonant multiple scattering: applications to localized plasmon excitations // JOSA A. - 2008. - Vol. 25. - P. 2549-2557.

14. Дьячков П.Н. Углеродные нанотрубки: строение, свойства, применение. - М.: БИНОМ: Лаборатория знаний, 2006.

15. Ziman J.M. Principles of the Theory of Solids. - London: Cambridge University Press, 1979.

16. Победря Б.Е. Лекции по тензорному анализу. - М.: Изд-во Моск. ун-та, 1986.

17. Wriedt T. Using the T-Matrix method for light scattering computations by non-axisymmetric particles: Superellipsoids and realistically shaped particles. Part // Part. Syst. Charact. - 2002. - Vol. 19. - P. 256-268.

Получено 09.09.2011