Точное и эффективное моделирование наномасштабного оптического захвата с использованием FETI-DP и тензора Максвелла

Аннотация

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

Введение

Плазмонные оптические пинцеты на основе поверхностных плазмонов (ПП) привлекают большое внимание и широко применяются для захвата наночастиц [1–6]. ПП — это резонансное явление, вызванное связью падающего света с определенной длиной волны и свободными электронами на границе раздела металлов и диэлектриков [7]. ПП позволяют оптическим пинцетам преодолеть дифракционный предел. Более того, резкое локальное усиление поля ПП может снизить требования к интенсивности падающего света [7, 8]. Однако ПП тесно связаны с материалом и размерами объектов, а также длиной волны падающего света, что требует большого числа экспериментов для определения оптимальных параметров ПП оптических пинцетов на практике. Исходя из этого, метод моделирования играет все более важную роль как вспомогательное средство для проектирования и оптимизации ПП оптических пинцетов [9]. В этих симуляциях требуется расчет оптической силы для предсказания стабильного захвата. Для регулярных объектов, таких как сферы, оптическая сила может быть аналитически выведена из обобщенной теории Лоренца-Ми [10, 11]. Однако для объектов со сложными конфигурациями необходимы численные методы, строго решающие уравнения Максвелла, для моделирования электромагнитных полей и последующей оптической силы и потенциала.

Эти численные методы можно в основном классифицировать на методы дифференциальных уравнений (МДУ) и методы интегральных уравнений (МИУ) [12–15]. По сравнению с МИУ, методы дифференциальных уравнений (МДУ) показывают превосходные способности в работе со сложными геометриями и компонентами. МДУ также имеют преимущество прямого расчета распределения ближнего поля, которое играет важную роль в анализе ПП. Как представитель МДУ, метод конечных разностей во временной области (FDTD) реализован во временной области, что позволяет легко получать полосовую информацию и переходные процессы [16, 17]. Однако FDTD требует точной дисперсионной модели для описания частотно-зависимых материаловых свойств в ПП, тогда как точность решения FDTD сильно зависит от точности аппроксимации этой дисперсионной модели [18]. Кроме того, FDTD опирается на структурированные сетки, что часто приводит к ошибкам «лесенки» для кривых поверхностей. Как другой представитель МДУ, метод конечных элементов (МКЭ) широко применяется, так как он легко обрабатывает дисперсивные материалы в частотной области и устраняет ошибку «лесенки» неструктурированной сеткой [19–22]. По сравнению с FDTD, МКЭ может напрямую использовать измеренные материальные параметры без какой-либо приближенной дисперсионной модели. Однако резкое локальное усиление поля в ПП требует тонких сеток в дискретизации МКЭ. Кроме того, объекты с большими размерами и множественные объекты драматически увеличивают число неизвестных. Эти факторы вызывают плохо обусловленные матричные системы и огромные вычислительные затраты, создающие большие вызовы для традиционного МКЭ для анализа ПП-усиленного оптического захвата.

В этой статье вводится эффективный метод двойной-первичной конечно-элементной разрыва и соединения (FETI-DP) для моделирования оптического захвата на наномасштабе. FETI-DP использует схему не перекрывающегося разложения области, которая делит исходную крупномасштабную сложную задачу на серию мелкомасштабных простых задач для их решения. Он обеспечивает условие передачи на интерфейсах подобластей для обеспечения непрерывности поля и вводит двойную переменную для редукции исходной трехмерной (3D) задачи к двумерной (2D) задаче через множитель Лагранжа. Первичные переменные в углах подобластей извлекаются для ускорения сходимости итерационного решения двойственной задачи [23–26]. Разработан подход разреженного представления низкого ранга для улучшения производительности FETI-DP. Он использует алгоритмы данных-разреженных представлений для улучшения эффективности решения подзадач и двойственной задачи [27, 28]. Предлагаемый метод обеспечивает полностью разделенные подобласти, позволяющие параллельное моделирование оптической силы для захвата наночастиц. Тензор напряжений Максвелла (МНМ), раскрывающий связь между электромагнитным полем и механическим импульсом, используется для оценки оптической силы [29]. На основе полученной оптической силы оптический потенциал может быть дополнительно рассчитан для анализа стабильного захвата. По сравнению с МИУ, предлагаемый метод более мощен в работе с составными материалами и решении ближнего поля для ПП-основанного оптического захвата. По сравнению с FDTD, предлагаемый метод может точно обрабатывать дисперсивные металлические материалы в системах ПП-основанного оптического захвата и устранять ошибку «лесенки» для объектов с кривыми границами. По сравнению с МКЭ, предлагаемый метод подходит для крупномасштабных вычислений оптического захвата. Несколько примеров проанализированы, и численные результаты демонстрируют точность и эффективность предлагаемого метода для предсказания и анализа оптического захвата на наномасштабе.

Методы

Формулировки FETI-DP

Для реализации FETI-DP исходная вычислительная область Ω сначала разрывается на серию неперекрывающихся подобластей Ωⁱ (i = 1, 2, 3⋯, N_s), как показано на рис. 1. В каждой подобласти Ωⁱ подобластная конечно-элементная система может быть получена из векторного волнового уравнения как

$$ \nabla \times {\mu}_r^{-1}\nabla \times {\mathbf{E}}^i-{k}_0^2{\varepsilon}_r{\mathbf{E}}^i={jk}_0{\eta}_0{\mathbf{J}}_{\mathrm{imp}}^i\kern1.08em \mathrm{in}\kern0.24em {\Omega}^i \quad (1) $$ $$ \hat{n}\times \nabla \times {\mathbf{E}}^i+{jk}_0\hat{n}\times \left(\hat{n}\times {\mathbf{E}}^i\right)=0\kern0.96em \mathrm{on}\kern0.24em {\Gamma}_{\mathrm{ABC}}^i \quad (2) $$

Схема разделения области с неперекрывающимися подобластями в методе FETI-DP. a Исходная область. b Разделенные подобласти и дискретизированные сетки.

где Eⁱ обозначает неизвестное электрическое поле, решаемое в Ωⁱ, k₀ и η₀ — волновое число и внутреннее сопротивление свободного пространства соответственно, а \( {\mathbf{J}}_i^{\mathrm{imp}} \) — наведенный ток. \( {\Gamma}_{\mathrm{ABC}}^i \) означает поглощающее граничное условие (ABC) для усечения бесконечной открытой области. Следует отметить, что k₀ должен быть заменен волновым сопротивлением в среде, если окружающая среда не является свободным пространством, что обычно для оптического захвата. На интерфейсе подобластей Γⁱ требуется предполагаемое граничное условие для генерации полной граничной задачи в Ωⁱ. Здесь накладывается условие передачи типа Робина с неизвестной вспомогательной переменной Λ_i как

$$ {\hat{n}}^i\times \left({\mu}_r^{-1}\nabla \times {\mathbf{E}}^i\right)+{\alpha}^i{\hat{n}}^i\times \left({\hat{n}}^i\times {\mathbf{E}}^i\right)={\boldsymbol{\Lambda}}^i\kern1.2em \mathrm{on}\kern0.36em {\Gamma}^i \quad (3) $$

где \( {\hat{n}}^i \) обозначает единичный внешний нормальный вектор на интерфейсе подобласти Γⁱ, а αⁱ — комплексный параметр, который часто можно выбрать как jk₀. Все подобласти затем дискретизируются тетраэдрическими элементами. В каждом элементе мы разворачиваем E с векторными базисными функциями N и неизвестными коэффициентами электрического поля E как

$$ \mathbf{E}=\sum \limits_{p=1}^s{E}_p{\mathbf{N}}_p \quad (4) $$

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

Применяя метод Галекина, матричное уравнение МКЭ в Ωⁱ относительно неизвестного коэффициента электрического поля Eⁱ может быть получено как

$$ \left(\begin{array}{cc}{\mathbf{K}}_{rr}^i& {\mathbf{K}}_{rc}^i\\ {}{\mathbf{K}}_{cr}^i& {\mathbf{K}}_{cc}^i\end{array}\right)\left(\begin{array}{l}{E}_r^i\\ {}{E}_c^i\end{array}\right)=\left(\begin{array}{l}{f}_r^i-{\mathbf{B}}_r^{i^T}{\lambda}^i\\ {}{f}_c^i\end{array}\right) \quad (5) $$

где индексы c и r отличают угловые степени свободы (DOF) и оставшиеся DOF, извлекая угловые DOF как первичную переменную для построения двойной-первичной (DP) схемы. Здесь K — системная матрица МКЭ, f — вектор возбуждения. B — булева матрица, извлекающая интерфейсные DOF подобласти. λ — двойственная переменная, генерируемая из развертывания Λ_i, также называемая множителем Лагранжа.

Затем смежные подобласти могут быть соединены путем обеспечения непрерывности касательных электрического и магнитного полей на интерфейсах. Мы собираем все интерфейсы подобластей и устраняем все внутренние неизвестные подобластей Eⁱ. Редуцированное глобальное интерфейсное уравнение относительно двойственной переменной λ может быть получено как

$$ \left[{\tilde{\mathbf{K}}}_{rr}+{\tilde{\mathbf{K}}}_{rc}{\tilde{\mathbf{K}}}_{cc}^{-1}{\tilde{\mathbf{K}}}_{cr}\right]\lambda ={\tilde{f}}_r-{\tilde{\mathbf{K}}}_{rc}{\tilde{\mathbf{K}}}_{cc}^{-1}{\tilde{f}}_c \quad (6) $$

Уравнение (6) может быть решено итерационными методами, такими как обобщенный метод минимального остатка (GMRES). \( {\tilde{\mathbf{K}}}_{cc} \) — глобальная угловая система, которая может ускорить итерационную сходимость в первичном пространстве. Подходящий предобуславливатель может улучшить скорость итерационной сходимости, такой как приближенная обратная и неполная LU-декомпозиция. После решения двойственной переменной λ, электрическое поле внутри каждой подобласти может быть независимо оценено по (5). Для построения глобальной матрицы \( {\tilde{\mathbf{K}}}_{rr} \), нужно построить подобластную численную функцию Грина \( {\mathbf{Z}}_{rr}^i \) в виде

$$ {\mathbf{Z}}_{rr}^i={\mathbf{B}}_r^i{\left({\mathbf{K}}_{rr}^i\right)}^{-1}{{\mathbf{B}}_r^i}^T \quad (7) $$

где включена обратная подобластная матрица МКЭ \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \). Кроме того, для матриц \( {\tilde{\mathbf{K}}}_{rc} \), \( {\tilde{\mathbf{K}}}_{cr} \), \( {\tilde{\mathbf{K}}}_{cc} \) и векторов \( {\tilde{f}}_r \), \( {\tilde{f}}_c \) требуется вычисление \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \). Конструкции \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) на этапе предобработки, а также их матрично-векторные произведения (MVP) на этапе итерационного решения являются вычислительно затратными. Хотя \( {\mathbf{K}}_{rr}^i \) разрежены, \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) плотны, что означает высокие вычислительные затраты. Далее вводится метод разреженного представления низкого ранга для ускорения вычислений \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \). Поскольку некоторые подматрицы в глобальной интерфейсной системе могут быть представлены в виде матриц низкого ранга, их вычисление можно выполнять с алгоритмами низкого ранга, что улучшает производительность FETI-DP. Видно, что FETI-DP обеспечивает независимые подобластные операции, поэтому он может использовать параллельные вычисления для повышения эффективности. Для эффективной параллельной схемы принцип разделения области — сделать число DOF в каждой подобласти максимально сбалансированным. Следовательно, размер подобластей должен относиться к плотности дискретизации сетки. Обычно малые подобласти применяются в тонко засетенных областях, а большие — в грубо засетенных.

Разреженное представление низкого ранга

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

$$ \mathbf{M}={\mathbf{XY}}^{\mathrm{T}}\kern0.72em \left(\mathbf{M}\in {\mathrm{\mathbb{C}}}^{m\times n},\mathbf{X}\in {\mathrm{\mathbb{C}}}^{m\times k},\mathbf{Y}\in {\mathrm{\mathbb{C}}}^{n\times k}\right) \quad (8) $$

где X и Y в полных матричных формах, а ранг k намного меньше m и n. Матрица \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) может быть представлена формами данных-разреженных матриц, так как она обладает матричным свойством интегрального оператора. Следовательно, при условии, что \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) обладает свойством низкого ранга в данной подобласти, она может быть эффективно вычислена и сохранена в формах данных-разреженных представлений с подходом разреженного представления низкого ранга, ускоряя MVP в итерационном решении.

Процессы подхода разреженного представления низкого ранга можно разделить на следующие шаги: (1) построить кластерное дерево путем подразделения множества базисных функций в каждой подобласти, (2) построить блочное кластерное дерево взаимодействием двух кластерных деревьев, (3) сгенерировать данные-разреженную форму \( {\mathbf{K}}_{rr}^i \) по условию допустимости, (4) выполнить алгоритмы низкорангового форматирования для получения данных-разреженного представления \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}}^{-1} \), и (5) войти в решение систем FETI-DP алгоритмом данных-разреженных представлений. Подходящий предобуславливатель может ускорить решение. Следует отметить, что данные-разреженная LU-факторизация \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}}={\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}}{\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}} \) используется для замены матричной инверсии \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}}^{-1} \). Техника вложенного сечения применяется для дальнейшего улучшения эффективности разреженного представления низкого ранга. Вложенное сечение использует разделители для получения больших внедиагональных нулевых подблоков, которые останутся нулевыми во время LU-факторизации, что значительно уменьшает заполнение.

Для генерации \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}} \), сначала строим кластерное дерево T_I рекурсивным подразделением множества реберных базисных функций подобласти I = {1,2,……N} с использованием ограничивающего параллелепипеда. С вложенным сечением кластер t внутри соответствующего ограничивающего параллелепипеда делится на трех преемников {s₁, s_sep, s₂}, где s₁ и s₂ — индексные множества двух неподключенных ограничивающих параллелепипедов, а s_sep — индексное множество разделителя. Рис. 2a показывает простой пример этого процесса. Затем блочное кластерное дерево T_{I × I} может быть построено взаимодействием двух кластерных деревьев T_I, как показано на рис. 2b, которое можно выбрать как кластерное дерево исходного реберного множества базисных функций и тестовых базисных функций в методе Галекина. Далее нужно ввести условие допустимости на основе вложенного сечения для различения полных блоков, блоков низкоранговой декомпозиции и внедиагональных нулевых блоков в T_{I × I} [23]. Таким образом, \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}} \) может быть произведен заполнением соответствующих блоков ненулевыми элементами \( {\mathbf{K}}_{rr}^i \). Наконец, данные-разреженная LU-факторизация \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}}={\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}}{\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}} \) может быть вычислена рекурсивно из

$$ {\mathbf{K}}_{rr}^i=\left[\begin{array}{ccc}{\mathbf{K}}_{11}& & {\mathbf{K}}_{13}\\ {}& {\mathbf{K}}_{22}& {\mathbf{K}}_{23}\\ {}{\mathbf{K}}_{31}& {\mathbf{K}}_{32}& {\mathbf{K}}_{33}\end{array}\right]=\left[\begin{array}{ccc}{\mathbf{L}}_{11}& & \\ {}& {\mathbf{L}}_{22}& \\ {}{\mathbf{L}}_{31}& {\mathbf{L}}_{32}& {\mathbf{L}}_{33}\end{array}\right]\left[\begin{array}{ccc}{\mathbf{U}}_{11}& & {\mathbf{U}}_{13}\\ {}& {\mathbf{U}}_{22}& {\mathbf{U}}_{23}\\ {}& & {\mathbf{U}}_{33}\end{array}\right] \quad (9) $$

Построение кластерного дерева и блочного кластерного дерева 4 уровней на основе вложенного сечения. a Построение кластерного дерева рекурсивным подразделением реберного множества базисных функций I = {1,2,…18}. b Построение блочного кластерного дерева, где белые блоки — нулевые матрицы, а зеленые блоки могут быть полными матрицами или матрицами низкоранговой декомпозиции.

где обычные полные матричные арифметические операции заменяются их данными-разреженными аналогами [28]. Адаптивная ошибка усечения ε_t используется для контроля точности низкоранговых аппроксимаций. Полученные LU-факторы \( {\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}} \) и \( {\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}} \) сохраняются и используются для построения \( {\mathbf{Z}}_{rr}^i \) по

$$ {\mathbf{Z}}_{rr}^i={\mathbf{B}}_r^i{\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}}^{-1}{\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}}^{-1}{\mathbf{B}}_r^i \quad (10) $$

где \( {\mathbf{B}}_r^i{\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}}^{-1} \) и \( {\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}}^{-1}{\mathbf{B}}_r^i \) могут быть вычислены данными-разреженными верхним и нижним треугольными решателями. \( {\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}} \), \( {\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}} \) и \( {\mathbf{Z}}_{rr}^i \) входят в расчет FETI-DP с данными-разреженными прямыми и обратными подстановками (FBS) и данными-разреженными MVP.

Оптическая сила и потенциал

Согласно теории электродинамики, оптическая сила может быть оценена тензором напряжений Максвелла (МНМ), раскрывающим связь между электромагнитным полем и механическим импульсом [29]. Получив распределение электромагнитного поля вокруг объекта, оптическую силу можно рассчитать интегрированием МНМ по замкнутой поверхности, окружающей объект. На основе полученного распределения электрического поля МНМ в любых координатах может быть построен как

$$ \overleftrightarrow{\mathbf{T}}=\frac{1}{2}\operatorname{Re}\left[\varepsilon {\mathbf{EE}}^{\ast }+\mu {\mathbf{HH}}^{\ast }-\frac{1}{2}\left(\varepsilon {\left|\mathbf{E}\right|}^2+\mu {\left|\mathbf{H}\right|}^2\right)\overleftrightarrow{\mathbf{I}}\right] \quad (11) $$

где верхний индекс звездочки обозначает комплексно-сопряженное электрическое или магнитное поле, ε и μ — диэлектрическая и магнитная проницаемость, а \( \overleftrightarrow{\mathbf{I}} \) — единичная матрица 3×3. Внешним произведением векторов тензорная форма \( \overleftrightarrow{\mathbf{T}} \) может быть записана как

$$ \overleftrightarrow{\mathrm{T}}=\left[\begin{array}{lll}{T}_{xx}& {T}_{xy}& {T}_{xz}\\ {}{T}_{yx}& {T}_{yy}& {T}_{yz}\\ {}{T}_{zx}& {T}_{zy}& {T}_{zz}\end{array}\right]=\left[\begin{array}{ccc}\varepsilon {E}_x{E}_x^{\ast }+\mu {H}_x{H}_x^{\ast }-\frac{\varepsilon {\left|\mathbf{E}\right|}^2+\mu {\left|\mathbf{H}\right|}^2}{2}& \varepsilon {E}_x{E}_y^{\ast }+\mu {H}_x{H}_y^{\ast }& \varepsilon {E}_x{E}_z^{\ast }+\mu {H}_x{H}_z^{\ast}\\ {}\varepsilon {E}_y{E}_x^{\ast }+\mu {H}_y{H}_x^{\ast }& \varepsilon {E}_y{E}_y^{\ast }+\mu {H}_y{H}_y^{\ast }-\frac{\varepsilon {\left|\mathbf{E}\right|}^2+\mu {\left|\mathbf{H}\right|}^2}{2}& \varepsilon {E}_y{E}_z^{\ast }+\mu {H}_y{H}_z^{\ast}\\ {}\varepsilon {E}_z{E}_x^{\ast }+\mu {H}_z{H}_x^{\ast }& \varepsilon {E}_z{E}_y^{\ast }+\mu {H}_z{H}_y^{\ast }& \varepsilon {E}_z{E}_z^{\ast }+\mu {H}_z{H}_z^{\ast }-\frac{\varepsilon {\left|\mathbf{E}\right|}^2+\mu {\left|\mathbf{H}\right|}^2}{2}\end{array}\right] \quad (12) $$

где индексы x, y, z обозначают компоненты в трех направлениях. Согласно развертыванию E, описанному в (4), элементы МНМ T_mn (m, n = x, y, z) могут быть преобразованы в развернутые формы в расчете FETI-DP как

$$ {\displaystyle \begin{array}{l}{T}_{mn}=\sum \limits_{p,q=1}^s{E}_p{E}_q\left\{\varepsilon {\left({\mathbf{N}}_p\right)}_m{\left({\mathbf{N}}_q^{\ast}\right)}_n-\frac{1}{\omega^2\mu }{\left(\nabla \times {\mathbf{N}}_p\right)}_m{\left(\nabla \times {\mathbf{N}}_q^{\ast}\right)}_n\right.\\ {}\kern1.75em \left.-\frac{1}{2}\left[\varepsilon \left({\mathbf{N}}_p\right)\left({\mathbf{N}}_q^{\ast}\right)-\frac{1}{\omega^2\mu}\left(\nabla \times {\mathbf{N}}_p\right)\left(\nabla \times {\mathbf{N}}_q^{\ast}\right)\right]\right\}\kern1.75em \mathrm{if}\ m=n.\end{array}} \quad (13) $$ $$ {T}_{mn}=\sum \limits_{p,q=1}^s{E}_p{E}_q\left\{\varepsilon {\left({\mathbf{N}}_p\right)}_m{\left({\mathbf{N}}_q^{\ast}\right)}_n-\frac{1}{\omega^2\mu }{\left(\nabla \times {\mathbf{N}}_p\right)}_m{\left(\nabla \times {\mathbf{N}}_q^{\ast}\right)}_n\right\}\kern1.25em \mathrm{if}\ m\ne n. \quad (14) $$

где ω — угловая частота; N и s описаны в Ур. (4).

Наконец, оптическая сила, действующая на объект, может быть рассчитана интегрированием МНМ по любой замкнутой поверхности, окружающей его, как

$$ \mathbf{F}={\oint}_S\left(\overleftrightarrow{\mathbf{T}}\cdot \hat{n}\right)\ dS. \quad (15) $$

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

Оптический потенциал — еще один привлекательный параметр, раскрывающий стабильность оптического захвата. На основе полученной оптической силы глубина оптического потенциала U в позиции r₀ может быть рассчитана как

$$ \mathbf{U}\left({r}_0\right)=-{\int}_{\infty}^{r_0}\mathbf{F}\left(\mathbf{r}\right)\cdot \mathbf{r}, \quad (16) $$

где индекс ∞ обозначает бесконечность, определенную как опорная точка с нулевым потенциалом. Значение U может быть представлено k_BT, где k_B — постоянная Больцмана 1,3806488 × 10⁻²³ Дж/К, T — окружающая температура. Обычно частица может преодолеть броуновское движение в растворе и быть стабильно захваченной, когда U > 1 k_BT выполняется. Иначе частица не может быть стабильно захвачена. Поскольку полная оптическая сила включает консервативную градиентную компоненту и неконсервативную рассеивающую компоненту, полная оптическая сила F из (15) неконсервативна [30, 31]. Однако, при условии, что движение наночастицы ограничено одним измерением, это дает однозначное определение оптического потенциала из (16), даже если полная оптическая сила неконсервативна.

Результаты и обсуждение

Представлены три примера для демонстрации эффективности предлагаемого метода. Поскольку благородные металлы обычно используются для возбуждения поверхностных плазмонов, мы выбираем репрезентативные золото и серебро для анализов. Первый пример вычисляет оптическую силу серебряной наночастицы для проверки точности предлагаемого метода. Второй и третий примеры моделируют и обсуждают оптический захват золотых наночастиц. Для всех примеров бесконечная область усечена ABC, а расстояния между ABC и объектами установлены как одна длина волны, что достаточно для достижения приемлемой точности. Все расчеты выполнены на рабочей станции Dell с процессорами Intel Xeon 3,6 ГГц.

Серебряная нанокапсула

Объект серебряной нанокапсулы сначала рассматривается для проверки точности и эффективности предлагаемого метода FETI-DP в предсказании оптической силы. Рис. 3a и b показывают его конфигурацию и размеры. Конститутивные параметры серебра — все измеренные значения из [32]. Для реализации схемы FETI-DP весь аналитический домен сначала делится на 24 подобласти. Плотные сетки требуются возле металлической поверхности для моделирования эффекта плазмонного локального усиления поля. Принимаются тетраэдрические элементы для дискретизации, что приводит к общему числу 6,9 × 10⁵ неизвестных, включая 4,1 × 10⁴ двойственных неизвестных и 313 угловых неизвестных. Падающее свет освещает вдоль направления +z, а направление поляризации электрического поля — −x.

Конфигурация структуры серебряной нанокапсулы. a 3D вид. b Фронтальный вид и размеры, где R = 30 нм и h = 60 нм.

Сначала меняем длину волны падающего света λ от 200 нм до 400 нм для моделирования оптических сил, действующих на нанокапсулу. Поскольку FETI-DP работает в частотной области, оптические силы рассчитываются в 15 точек выборки частоты. Рис. 4 показывает рассчитанную кривую оптических сил, действующих на серебряную нанокапсулу. Для указания точности FETI-DP результаты оптической силы FETI-DP сравниваются с результатами коммерческого ПО Lumerical FDTD Solutions [33], и наблюдается хорошее согласование.

Результаты оптических сил, действующих на серебряную нанокапсулу, зависящие от длины волны λ падающего света, включая результаты FETI-DP и коммерческого ПО FDTD Solutions.

Затем производительность FETI-DP тестируется для разных чисел подобластей. Увеличиваем число подобластей от 4 до 24, сохраняя плотность дискретизации. Каждому процессору назначаем одну подобласть. Таблица 1 сообщает время, затраченное на построение глобального интерфейсного Ур. (6) и общее время решения. Видно, что FETI-DP может полностью использовать параллельные вычислительные ресурсы и значительно улучшить эффективность решения. Кроме того, точность FETI-DP с увеличением числа подобластей также проверяется и сообщается в Таблице 1. Здесь точность определяется 2-нормой относительной ошибки оптической силы как δ_OF = ‖OF_i − OF_ref‖/‖OF_ref‖, где OF_i — оптическая сила с i подобластями, OF_ref — эталонная оптическая сила с двумя подобластями. Видно, что точность остается почти постоянной с увеличением числа подобластей.

Димер золотых наносфер

Второй пример анализирует оптический захват золотой наносферы с использованием димера золотых наносфер. Плазмонные эффекты в зазоре димера могут эффективно усилить оптическую силу для захвата наночастицы. Рис. 5a и b дают конфигурацию и размеры этой системы. Конститутивные параметры золота — все измеренные значения из [32]. Окружающая среда — вода с относительным показателем преломления n = 1,33. Падающий свет — плоская волна с мощностью 10 мВт/мкм², поляризация электрического поля +x, направление падения −z. Оптическая сила, действующая на объект наносферу, рассчитывается методом FETI-DP. Для реализации FETI-DP весь вычислительный домен делится на 32 подобласти и дискретизируется тетраэдрическими сетками, что приводит к 3,5 × 10⁶ неизвестных, включая 1,6 × 10⁵ двойственных неизвестных и 1738 угловых неизвестных.

Конфигурация системы оптического захвата димера золотых наносфер в воде. a 3D вид. b Фронтальный вид и размеры, где R = 25 нм, r = 5 нм, g = 2 нм.

Сначала тестируем параллельную производительность предлагаемого FETI-DP с использованием различных чисел процессоров. Таблица 2 сообщает время решения для Ур. (6) и общее время решения. Кроме того, ускорения для параллельных вычислений также даны в Таблице 2. Здесь ускорение определяется как Speed up = T₁/T_Np, где T_Np обозначает общее настенное время с N_p процессорами. Видно, что FETI-DP значительно улучшает эффективность решения и демонстрирует хорошее параллельное ускорение. Для этого большого числа неизвестных общее использование памяти всех процессоров составляет только 57,2 ГБ.

Затем проверяется эффективность подхода разреженного представления низкого ранга. С разреженным представлением низкого ранга подобластная матрица может быть факторизована алгоритмом данных-разреженных представлений и сохранена как данные-разреженные матрицы. Время построения и использование памяти составляют только 18 с и 0,5 ГБ, тогда как они 67 с и 1,7 ГБ при обычном матричном алгоритме. Видно, что мы получаем 72% экономию времени и 70% сжатие памяти. Что касается использования памяти, последующие MVP также получают 70% экономию времени.

Далее FETI-DP тестируется для расчета оптической силы с длиной волны λ, варьирующей от 277 нм до 818 нм. На практике анализы оптической силы под падающим светом разных длин волны часто необходимы для поиска плазмонной резонансной длины волны, где происходит резкое усиление поля и может быть получена максимальная оптическая сила. Рассматриваются два случая с наносферой, расположенной в (0, 0, 20 нм) и (0, 0, −20 нм). Рис. 6a и b строят рассчитанные оптические силы, действующие на наносферу для разных λ. Видно, что максимальная оптическая сила происходит при λ = 472 нм, что является плазмонной резонансной длиной волны. Оптическая сила на этой резонансной длине волны усилена почти в 40 раз по сравнению с нерезонансной длиной волны. Более того, оптическая сила всегда направлена в зазор димера, как показано на рис. 6, где интенсивность электрического поля максимальна. Это также направление градиентной силы для захвата объекта. Рис. 7a и b показывают рассчитанные распределения усиления электрического поля на нерезонансной длине волны λ = 300 нм и резонансной длине волны λ = 472 нм соответственно. Видно, что интенсивность электрического поля увеличена почти в 250 раз благодаря плазмонному резонансному эффекту.

Рассчитанные результаты оптических сил, действующих на наносферу в системе димера золотых наносфер, зависящие от длины волны λ падающего света. a Объект наносфера находится в (0, 0, 20 нм). b Объект наносфера находится в (0, 0, −20 нм).

Распределения усиления электрического поля на плоскости xoz для системы димера золотых наносфер. a λ = 300 нм (нерезонансная длина волны). b λ = 472 нм (резонансная длина волны).

Кроме того, оптическая сила и оптический потенциал рассчитываются с наносферой, движущейся от (0, 0, −30 нм) до (0, 0, −17 нм) вдоль оси z. Поскольку наиболее типичное и интересное поведение сил и потенциалов захвата — те, действующие вдоль z-направления, мы здесь рассматриваем осевой захват потенциала интегрированием вдоль z-оси. Поскольку движение наночастицы ограничено одним измерением, определение оптического потенциала из (16) однозначно, даже если полная оптическая сила из (15) неконсервативна. Как показано на рис. 8a, b, с движением наносферы к зазору димера оптическая сила и глубина оптического потенциала очевидно увеличиваются. В позиции (0, 0, −17 нм) получена глубина оптического потенциала 4,6 k_BT, достаточная для преодоления броуновского движения в воде для достижения стабильного оптического захвата.

Оптические силы и оптические потенциалы, действующие на наносферу в системе димера золотых наносфер, когда наносфера движется от (0, 0, −30 нм) до (0, 0, −17 нм). a Оптические силы. b Оптические потенциалы.

Наконец, тестируем эффекты диэлектрического подложки для этого примера. Оптические силы рассчитывают с и без подложки соответственно. Для обоих случаев наносфера находится в (0, 0, −20 нм) и длина волны падения выбрана как резонансная длина волны. Для случая без подложки рассчитанный результат оптической силы |F₀| = 0,769 пН. Для случая с подложкой димер золотых наносфер помещен на диэлектрическое подложку толщиной 60 нм и относительной диэлектрической проницаемостью ε_r = 2,25. Рассчитанный результат оптической силы |F₁| = 0,761 пН. Относительная ошибка между этими двумя результатами оптических сил составляет около 1,0 × 10⁻², определенная как |F₁ − F₀|/|F₀|.

Димер золотых усеченных конусов

Третий пример касается оптического захвата золотой наносферы с использованием димера золотых усеченных конусов. Рис. 9 дает конфигурацию и размеры этой системы. Конститутивные параметры золота взяты из [32]. Диэлектрическое подложка имеет относительную диэлектрическую проницаемость ε_r = 2,25. Окружающая среда — вода с относительным показателем преломления n = 1,33. Падающий свет — плоская волна с мощностью 10 мВт/мкм², поляризация электрического поля +x, направление падения −z. Весь вычислительный домен делится на 32 подобласти и дискретизируется тетраэдрическими сетками, что приводит к 3,1 × 10⁶ неизвестных, включая 1,3 × 10⁵ двойственных неизвестных и 1227 угловых неизвестных.

Конфигурация системы оптического захвата димера золотых усеченных конусов на диэлектрическом подложке в воде. a 3D вид. b Фронтальный вид и размеры, где UR = 20 нм, LR = 30 нм, h = 35 нм, g = 2 нм.

Сначала анализируем оптические силы, меняя λ от 277 нм до 818 нм. Рис. 10 строит рассчитанные оптические силы, действующие на наносферу для разных λ. Наносфера находится в (0, 0, 35 нм). Видно, что максимальная оптическая сила происходит при λ = 464 нм, что является плазмонной резонансной длиной волны, и оптическая сила здесь усилена почти в 30 раз при нерезонансной длине волны. Более того, полная оптическая сила всегда направлена в −z, как показано на рис. 10, что является направлением градиентной силы. Это подтверждает, что градиентная сила больше рассеивающей силы, что является одним из условий стабильного захвата наносферы. Рис. 11a и b представляют рассчитанные распределения электрического поля на нерезонансной длине волны λ = 300 нм и резонансной длине волны λ = 464 нм соответственно. Видно, что интенсивность электрического поля увеличена почти в 500 раз благодаря локализованному поверхностному плазмонному резонансу.

Рассчитанные результаты оптических сил, действующих на наносферу в системе димера золотых усеченных конусов, зависящие от λ. Наносфера находится в (0, 0, 35 нм).

Распределения усиления электрического поля на плоскости xoz для системы димера золотых усеченных конусов. a λ = 300 нм (нерезонансная длина волны). b λ = 464 нм (резонансная длина волны).

Затем положение наносферы меняется на 0, 5 и 35 нм для наблюдения оптической силы. Рис. 12 дает рассчитанные оптические силы, действующие на наносферу, где наблюдается очевидная y-компонента оптической силы, в то время как существует большая z-компонента. Полная оптическая сила все равно указывает в позицию с сильнейшим электрическим полем для захвата наносферы.

Рассчитанные результаты оптических сил, действующих на наносферу в системе димера золотых усеченных конусов, варьирующие λ. Наносфера находится в (0, 5 нм, 35 нм).

Далее анализируем оптический потенциал с наносферой, движущейся от (0, 0, 50 нм) до (0, 0, 20 нм) вдоль оси z. Здесь рассматриваем осевой захват потенциала вдоль z-направления, ограничивающее движение наночастицы одним измерением и ведущее к однозначному определению оптического потенциала. И оптическая сила, и потенциал рассчитываются. Как видно из рис. 13a, b, с движением наносферы к зазору димера оптическая сила и глубина оптического потенциала очевидно увеличиваются. В (0, 0, 20 нм) получена глубина оптического потенциала 3,8 k_BT, достаточная для преодоления броуновского движения в воде для достижения стабильного оптического захвата.

Оптические силы и оптические потенциалы, действующие на наносферу в системе димера золотых усеченных конусов, когда наносфера движется от (0, 0, 50 нм) до (0, 0, 20 нм). a Оптические силы. b Оптические потенциалы.

Наконец, тестируем вычислительные затраты FETI-DP, меняя число неизвестных от 1,0 миллиона до 3,2 миллиона на основе разных размеров сетки. На практике тесты под разной плотностью сетки обычно необходимы для удовлетворения разных требований точности. Такая крупномасштабная сложная задача создает большие вызовы для традиционных численных методов. Однако FETI-DP может легко обработать эту проблему. Тридцать два процессора используются для симуляции FETI-DP, каждый процессор работает с подобластью. Таблица 3 сообщает вычислительные затраты FETI-DP. Видно, что FETI-DP демонстрирует высокую эффективность моделирования и низкие требования к памяти.

Заключение

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

Доступность данных и материалов

Все данные, сгенерированные или проанализированные в этом исследовании, включены в эту статью.

Аббревиатуры

ABC:

Поглощающее граничное условие

DOF:

Степени свободы

FDTD:

Метод конечных разностей во временной области

FEM:

Метод конечных элементов

FETI-DP:

Двойной-первичный метод конечно-элементного разрыва и соединения

MST:

Тензор напряжений Максвелла

MVP:

Матрично-векторное произведение