Лицензия MIT
© Алексей Александрович Щербаков, 2024
Теория Ми позволяет эффективно рассчитывать параметры рассеяния на сферических частицах, а её обобщения - на многослойных сферических частицах. Если форма рассеивателя отличается от сферической, расчет существенно усложняется. Для несферических рассеивателей был предложен ряд специальных методов, включающих приближение дискретных диполей, метод продолженных граничных условий, метод дискретных источников, метод инвариатного вложения, а также их различные модификации.
В данной лекции мы рассмотрим метод дискретных диполей (Discrete Dipole Approximation), часто применяющийся при решении научных и прикладных задач и имеющий ряд эффективных реализаций в форме открытого кода. Метод основан на дискретизации объема и в разных источниках можно встретить его упоминания как метода моментов (для разбиения объема) (Method of Moments), метода тензорных функций Грина (Green's dyadic formalism), или просто метода объемного интегрального уравнения (Volume Integral Equation).
Мы будем рассматривать рассеивающие частицы, занимающие пространственную область $V_s$ и расположенные в однородной изотропной среде, так что диэлектрическая проницаемость во всем пространстве задана функцией $$ \varepsilon(\boldsymbol{r}) = \begin{cases} \varepsilon_m, & \boldsymbol{r}\notin V_s \\ \varepsilon(\boldsymbol{r}), & \boldsymbol{r}\in V_s \end{cases} $$ Для простоты будем считать магнитную проницаемость равной проницаемости вакуума $\mu_0$ всюду. Обобщение метода на случай анизотропных, магнитных и бианизотропных сред не представляет принципиальных сложностей.
Приближение дискретных диполей основано на объемном интегральном уравнении с тензорной функцией Грина свободного пространства, являющегося решением векторного уравнения Гельмгольца с волновым числом в свободном пространстве, окружающем частицу, $k_m = \omega\sqrt{\varepsilon_m\mu_0}$. Рассмотрим частицу как изменение диэлектрической проницаемости однородного пространства $\Delta\varepsilon=\varepsilon(\boldsymbol{r})-\varepsilon_m$. Тогда наличие частицы можно трактовать как эффективный ток поляризации \begin{equation} \boldsymbol{J}=-i\omega\Delta\varepsilon\boldsymbol{E}. \end{equation} и пользоваться тензорной функцией Грина свободного пространства. Запишем уравнение для интеграла с вмысле главного значения, как было продемонстрировано во вводной лекции по функциям Грина \begin{equation}\tag{1} \boldsymbol{E}\left(\boldsymbol{r}\right) = \boldsymbol{E}^{inc}\left(\boldsymbol{r}\right) + k_{m}^{2}\lim_{V_0\rightarrow 0} \intop_{V_s\setminus V_0} \mathcal{G}_{0} \left(\boldsymbol{r}-\boldsymbol{r}'\right) \dfrac{\Delta\varepsilon\left(\boldsymbol{r}'\right)}{\varepsilon_{m}}\boldsymbol{E}\left(\boldsymbol{r}'\right)d^{3}\boldsymbol{r}' - \mathcal{L}(\boldsymbol{r}, \partial V_0) \dfrac{\Delta\varepsilon\left(\boldsymbol{r}\right)}{\varepsilon_{m}}\boldsymbol{E}\left(\boldsymbol{r}\right) \end{equation} Здесь $\boldsymbol{E}^{inc}$ - известное падающее поле, а тензорная функция Грина в явном виде \begin{equation}\tag{2} \mathcal{G}_0(\boldsymbol{R}) = \dfrac{e^{ikR}}{4\pi R} \left[ \mathbb{I}\left(1+\dfrac{i}{kR}-\dfrac{1}{k^{2}R^{2}}\right)+\hat{\boldsymbol{e}}_r \hat{\boldsymbol{e}}_r^T \left(\dfrac{3}{k^{2}R^{2}} - \dfrac{3i}{kR} - 1\right)\right] \end{equation} где $\boldsymbol{R} = \boldsymbol{r}-\boldsymbol{r}'$, $R = |\boldsymbol{R}|$, $\hat{\boldsymbol{e}}_r = \boldsymbol{R}/R$. Объем $V_0$ окружает точку $\boldsymbol{r}$. Тензор $\mathcal{L}$ иногда называется тензором деполяризации и имеет следующий явный вид, зависяций о формы поверхности исключаемого объема $V_0$: \begin{equation}\tag{2} \mathcal{L}\left(\boldsymbol{r},\partial V_0\right) = \oint_{\partial V_0}d^{2} \boldsymbol{r}' \hat{\boldsymbol{n}}' \dfrac{\left(\boldsymbol{r}'-\boldsymbol{r}\right)^{T}}{4\pi\left|\boldsymbol{r}'-\boldsymbol{r}\right|^{3}} \end{equation}
Для численного решения уравнения Липпмана-Швингера разобъем объем частицы на электрически малые непересекающиеся объемные ячейки $V_m$, как правило, кубической формы, с одинаковым объёмом $\Delta V = a^3$, так что $V_s=\bigcup_{k=1}^{N}V_{k}$, и центры ячеек оказываются расположены на регулярной сетке $\boldsymbol{r}_k = ak_x\hat{\boldsymbol{e}}_x + ak_y\hat{\boldsymbol{e}}_y + ak_z\hat{\boldsymbol{e}}_z$, $k_{x,y,z}\in\mathbb{Z}$, то есть под $k$ подразумевается трехмерный индекс для нумерации вдоль трёх пространственных направлений. Тогда, считая поле постоянным в объёме каждой ячейки, (что соответствует методу Галеркина с дельта-функциями $\delta(\boldsymbol{r}-\boldsymbol{r}_k)$, взятых в качестве базовых и тестовых), получаем следующее уравнение на поля в ячейках: \begin{equation}\tag{3} \boldsymbol{E}\left(\boldsymbol{r}_k\right) = \boldsymbol{E}^{inc}\left(\boldsymbol{r}_k\right) + k_{m}^{2} \Delta V \sum_{l=1,l\neq k}^{N} \mathcal{G}_{kl} \dfrac{\Delta\varepsilon_l}{\varepsilon_{m}} \boldsymbol{E}_l + (\mathcal{M}_k-\mathcal{L}_k) \dfrac{\Delta\varepsilon_k}{\varepsilon_{m}} \boldsymbol{E}_k \end{equation} Здесь введены обозначения $\boldsymbol{E}_k = \boldsymbol{E}\left(\boldsymbol{r}_k\right)$, $\Delta\varepsilon_k = \Delta\varepsilon(\boldsymbol{r}_k)$, и $\mathcal{G}_{kl} = \mathcal{G}_{0} \left(\boldsymbol{r}_k-\boldsymbol{r}_l\right)$. Тензор \begin{equation}\tag{4} \mathcal{M}_k\left(\boldsymbol{r},\Delta V\right) = k_m^2 \int_{\Delta V}\left[\mathcal{G}_{0}\left(\boldsymbol{r}_k-\boldsymbol{r}'\right)-\mathcal{G}_{s}\left(\boldsymbol{r}_{k}-\boldsymbol{r}'\right)\right] d^{3}\boldsymbol{r}' \end{equation} можно найти численно, но в расчетах, когда размер $a\ll\lambda$, им зачастую вообще пренебрегают. Вычитание статической части тензора Грина $\mathcal{G}_{s} = \lim_{kR\rightarrow 0} \mathcal{G}_{0}$, которая в явном виде записывается как $\mathcal{G}_{s} = -(4\pi k^2R^3)^{-1}(\mathbb{I}-3\hat{\boldsymbol{e}}_R\hat{\boldsymbol{e}}_R^T)$, делает подынтегральное выражение слабо сингулярным.
Уравнения (3) задают систему $3N$ самосогласованных линейных уравнений на неизвестные поля в каждой ячейке. Систему можно записать в матрично-векторном виде как \begin{equation}\tag{5} \left(\mathbb{I} - \mathrm{A}\right) \boldsymbol{u} = \boldsymbol{u}^{inc} \end{equation} где правая часть $\boldsymbol{u}^{inc} = \{ E^{inc}_{\alpha,k} \}_{k=1}^N$, $\alpha=x,y,z$, неизвестный вектор $\boldsymbol{u} = \{ E_{\alpha,k} \}_{k=1}^N$, а элементы матрицы $A$ задаются уравнением (3): \begin{equation}\tag{6} \mathrm{A}_{kl} = k_m^2 \Delta V \dfrac{\Delta\varepsilon_{l}}{\varepsilon_m}\mathcal{G}_{kl} + \delta_{kl}\dfrac{\Delta\varepsilon_k}{\varepsilon_m} \left(\mathcal{M}_k - \mathcal{L}_k\right) \end{equation}
В представленной формулировке метод представляет собой метод моментов. После решения системы (5) поле в произвольной точке в ближней зоне рассеянного поля можно найти с помощью прямого суммирования: \begin{equation}\tag{7} \boldsymbol{E}\left(\boldsymbol{r}\right) = \boldsymbol{E}^{inc}\left(\boldsymbol{r}\right) + \Delta V k_{m}^{2} \sum_{k=1}^{N} \mathcal{G}_{0}\left(\boldsymbol{r},\boldsymbol{r}_{k}\right) \dfrac{\Delta\varepsilon_k}{\varepsilon_m} \boldsymbol{E}_k,\thinspace\boldsymbol{r}\notin V_s \end{equation} Для нахождения амплитуды поля в дальней зоне в направлении, заданном вектором $ \hat{\boldsymbol{e}}_r$, следует воспользоваться асимптотикой функции Грина \begin{equation*} \mathcal{G}_{0}(\boldsymbol{r}, \boldsymbol{r}') \sim \left(\mathbb{I} - \hat{\boldsymbol{e}}_r \hat{\boldsymbol{e}}_r^T\right)\dfrac{e^{ikr}}{4\pi r}e^{-ik\hat{\boldsymbol{e}}_r\cdot\boldsymbol{r}'},\;r\gg r' \end{equation*} что даёт представление рассеянного поля \begin{equation}\tag{8} \boldsymbol{E}^{sca}\left(\hat{\boldsymbol{e}}_r\right) \sim \dfrac{\exp\left(ik_{b}r\right)}{r}\boldsymbol{F}^{sca}\left(\hat{\boldsymbol{e}}_r\right) \end{equation} с амплитудой рассеяния \begin{equation}\tag{9} \boldsymbol{F}^{sca}\left(\hat{\boldsymbol{e}}_r\right) = -k_m^{2} \dfrac{\Delta V}{4\pi} \sum_{k} \dfrac{\Delta\varepsilon_k}{\varepsilon_m}\exp\left(-ik_m\hat{\boldsymbol{e}}_r\cdot\boldsymbol{r}_k\right) \hat{\boldsymbol{e}}_r\times\hat{\boldsymbol{e}}_r\times\boldsymbol{E}_k \end{equation}
В рамках формулировки метода дискретных диполей вместо того, чтобы решать уравнение на неизвестные поля, можно выделить, так называемое, возбуждающее поле в каждой ячейке, вычитая диагональный член, соответствующий самодействию в формуле (3): \begin{equation}\tag{10} \boldsymbol{E}_k^{exc} = \boldsymbol{E}_k - \boldsymbol{E}_k^{self} = \boldsymbol{E}_k-\dfrac{\Delta\varepsilon_k}{\varepsilon_m}\left(k_m^{2}\mathcal{M}_k-\mathcal{L}_k\right)\boldsymbol{E}_k \end{equation} Тогда самосогласованная система на это поле запишется как \begin{equation}\tag{11} \boldsymbol{E}_k^{inc} = \sum_l \left(\delta_{kl}+\mathcal{G}_{kl} \boldsymbol{\alpha}_l\right) \boldsymbol{E}_l^{exc} \end{equation} где введен тензор поляризуемости \begin{equation}\tag{12} \boldsymbol{\alpha}_k = \dfrac{\Delta\varepsilon_k}{\varepsilon_m}\Delta V\left[\mathbb{I}-\dfrac{\Delta\varepsilon_k}{\varepsilon_m}\left(k_m^2\mathcal{M}_k-\mathcal{L}_k\right)\right]^{-1} \end{equation} связывающий вектор поляризации и возбуждающее поле в каждой ячейке: \begin{equation}\tag{13} \boldsymbol{P}_{k} = \boldsymbol{\alpha}_{k} \boldsymbol{E}_{k}^{exc} \end{equation} Заметим, что обращение в (13) относится к диагональной матрице, и, поэтому, не представляет вычислительной сложности. Итоговое уравнение можно сформулировать как уравнение на неизвестные поляризации: \begin{equation}\tag{14} \sum_{l}\left( \delta_{kl}\boldsymbol{\alpha}_{k}^{-1}+\mathcal{G}_{kl} \right)\boldsymbol{P}_{l} = \boldsymbol{E}_{k}^{inc} \end{equation} В такой формулировке подход к расчету полей назвается приближением дискретных диполей. Название связано с тем, что простейшая формулировка метода может быть получена исходя из предположения, что отклик каждой малой объемной ячейки является дипольным, и требуется рассчитать взаимодействие всех заданных диполей друг с другом и с внешним полем.
Сечения можно найти с помощью аплитуды рассеяния (9). Дифференциальная мощность рассеяния \begin{equation}\tag{15} \frac{dP^{sca}}{d\Omega} = \frac{1}{2Z_m} |\boldsymbol{F}^{sca}(\hat{\boldsymbol{e}}_r)|^2 \end{equation} Полное сечение рассеяния \begin{equation}\tag{16} C_{sca} = \frac{1}{k_m^2} \oint_{4\pi} |\boldsymbol{F}^{sca}(\hat{\boldsymbol{e}}_r)|^2 d\Omega = \frac{k_m^3}{E_{inc}^2} \sum_{k,l} \boldsymbol{P}_k^* \cdot \Im m \{ \mathcal{G}_{kl} \} \cdot \boldsymbol{P}_l \end{equation} где $E_{inc}$ - амплитуда поля падающей плоской волны. Сечение поглощения \begin{equation}\tag{17} C_{abs} = -\frac{k_m}{E_{inc}^2} \sum_{k} \Im m \{ \boldsymbol{E}_k^* \cdot \boldsymbol{P}_k^* \} \end{equation} и сечение экстинкции \begin{equation}\tag{18} C_{ext} = -\frac{k_m}{E_{inc}^2} \sum_{k} \Im m \{ \boldsymbol{E}_{inc,k}^* \cdot \boldsymbol{P}_k^* \} \end{equation} где $ \boldsymbol{E}_{inc,k}$ - векторная комплексная амплитуда поля падающей плоской волны в $k$-й ячейке.