Simulación precisa y eficiente de atrapamiento óptico nanoescala usando FETI-DP y tensor de Maxwell

Resumen

La simulación numérica juega un papel importante para la predicción de atrapamiento óptico basado en pinzas nano-ópticas plasmonicas. Sin embargo, estructuras complicadas y mejora drástica de campo local de efectos plasmonicos traen grandes desafíos a métodos numéricos tradicionales. En este artículo, se propone un método de simulación numérica preciso y eficiente basado en un método dual-primal de desgarro e interconexión de elementos finitos (FETI-DP) y tensor de estrés de Maxwell, para calcular la fuerza y potencial óptico para atrapar nanopartículas. Se introduce un enfoque de esparsificación de bajo rango paraさらに mejorar el rendimiento de simulación FETI-DP. El método propuesto puede descomponer un problema complejo y de gran escala en problemas simples y de pequeña escala usando división de dominio no superpuesta y discretización de malla flexible, lo que exhibe alta eficiencia y paralelizabilidad. Resultados numéricos muestran la efectividad del método propuesto para la predicción y análisis de atrapamiento óptico a nanoescala.

Introducción

Las pinzas ópticas plasmonicas basadas en plasmones de superficie (SPs) atraen mucha atención y se han aplicado ampliamente para capturar nanopartículas [1-6]. SP es un fenómeno de resonancia causado por el acoplamiento de luz incidente con una longitud de onda específica y electrones libres en la interfaz de metales y dieléctricos [7]. Los SPs permiten a las pinzas ópticas romper el límite de difracción. Además, la mejora drástica de campo local de los SPs puede reducir la demanda de intensidad de luz incidente [7, 8]. Sin embargo, los SPs están estrechamente relacionados con el material y dimensiones de objetos así como la longitud de onda de luz incidente, lo que requiere un gran número de experimentos para determinar los parámetros óptimos de pinzas ópticas SP en la práctica. Basado en esto, el método de simulación juega un papel cada vez más importante como medio auxiliar para el diseño y optimización de pinzas ópticas SP [9]. En estas simulaciones, el cálculo de fuerza óptica es requerido para predecir un atrapamiento estable. Para objetos regulares como esferas, la fuerza óptica puede derivarse analíticamente de la teoría generalizada de Lorenz-Mie [10, 11]. Sin embargo, para objetos con configuraciones complicadas, métodos numéricos que resuelven rigurosamente las ecuaciones de Maxwell son necesarios para modelar los campos electromagnéticos y la subsecuente fuerza y potencial óptico.

Estos métodos numéricos pueden categorizarse principalmente en métodos de ecuaciones diferenciales (DEMs) y métodos de ecuaciones integrales (IEMs) [12-15]. Comparado con IEMs, los métodos de ecuaciones diferenciales (DEMs) muestran habilidades superiores para lidiar con geometrías y componentes complicados. Los DEMs también tienen la ventaja de un cálculo directo de distribución de campo cercano, que juega un papel importante en el análisis de SPs. Como un DEM representativo, el método de diferencias finitas en dominio temporal (FDTD) se implementa en el dominio temporal, que puede obtener fácilmente información de banda ancha y respuestas transitorias [16, 17]. Sin embargo, el FDTD demanda un modelo dispersivo preciso para describir las propiedades de material dependientes de frecuencia en SPs, mientras que la precisión de la solución FDTD depende altamente de la precisión de aproximación de este modelo dispersivo [18]. Además, el FDTD depende de mallas estructuradas, que a menudo llevan a error de escalera para superficies curvas. Como otro DEM representativo, el método de elementos finitos (FEM) ha sido ampliamente adoptado ya que puede manejar fácilmente material dispersivo en el dominio de frecuencia y eliminar el error de escalera por malla no estructurada [19-22]. Comparado con FDTD, el FEM puede adoptar directamente parámetros de material medidos sin ningún modelo dispersivo aproximado. Sin embargo, las mejoras drásticas de campo local en los SPs requieren mallas finas en la discretización FEM. Además, objetos con grandes dimensiones y múltiples objetos aumentarán dramáticamente el número de incógnitas. Estos factores causarán sistemas de matrices mal condicionados y enormes consumos computacionales, trayendo grandes desafíos al FEM tradicional para el análisis de atrapamiento óptico mejorado por SP.

En este artículo, se introduce un eficiente método dual-primal de desgarro e interconexión de elementos finitos (FETI-DP) para simular el atrapamiento óptico a nanoescala. El FETI-DP adopta un esquema de descomposición de dominio no superpuesto, que divide un problema complejo original de gran escala en una serie de problemas simples de pequeña escala para conquistarlos. Impone una condición de transmisión en las interfaces de subdominio para asegurar la continuidad del campo, e introduce una variable dual para reducir el problema tridimensional (3D) original a un problema bidimensional (2D) por multiplicador de Lagrange. Variables primales en las esquinas de subdominio se extraen para acelerar la tasa de convergencia de la solución iterativa del problema dual [23-26]. Se desarrolla un enfoque de esparsificación de bajo rango para mejorar el rendimiento del FETI-DP. Usa algoritmos de datos dispersos para mejorar la eficiencia para resolver los problemas de subdominio y el problema dual [27, 28]. El método propuesto proporciona subdominios completamente desacoplados, que permiten la simulación paralela de fuerza óptica para atrapar nanopartículas. El tensor de estrés de Maxwell (MST) que revela la relación entre el campo electromagnético y momento mecánico se adopta para evaluar la fuerza óptica [29]. Basado en la fuerza óptica obtenida, el potencial óptico puede calcularseさらに para el análisis de un atrapamiento estable. Comparado con IEMs, el método propuesto es más poderoso para lidiar con materiales compuestos y resolver el campo cercano para atrapamiento óptico basado en SP. Comparado con FDTD, el método propuesto puede manejar con precisión material metálico dispersivo en sistemas de atrapamiento óptico basados en SP y eliminar el error de escalera para objetos con frontera curva. Comparado con FEM, el método propuesto es adecuado para computación de gran escala de atrapamiento óptico. Varios ejemplos se analizan y resultados numéricos demuestran la precisión y eficiencia del método propuesto para la predicción y análisis de atrapamiento óptico a nanoescala.

Métodos

Formulaciones FETI-DP

Para la implementación FETI-DP, el dominio computacional original Ω se desgarra primero en una serie de subdominios no superpuestos Ωⁱ (i = 1, 2, 3⋯, N_s), como se muestra en Fig. 1. En cada subdominio Ωⁱ, un sistema de elementos finitos de subdominio puede derivarse de la ecuación de onda vectorial como

$$ \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 $$ (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 $$ (2)

Esquema de división de dominio con subdominios no superpuestos en el método FETI-DP. a Dominio original. b Subdominios divididos y mallas discretizadas

donde Eⁱ denota el campo eléctrico desconocido a resolver en \( {\Omega}^i \), k₀ y η₀ son el número de onda en espacio libre e impedancia intrínseca, respectivamente, y \( {\mathbf{J}}_i^{\mathrm{imp}} \) es la corriente impresa. \( {\Gamma}_{\mathrm{ABC}}^i \) significa la condición de frontera absorbente (ABC) para truncar la región abierta infinita. Debe notarse que k₀ debe ser reemplazado por la impedancia de onda en medio si el medio circundante no es espacio libre, lo cual es común para el atrapamiento óptico. En la interfaz de subdominio Γⁱ, se requiere una condición de frontera asumida para generar un problema de valor de frontera completo en Ωⁱ. Aquí, se impone una condición de transmisión tipo Robin con una variable auxiliar desconocida Λ_i como

$$ {\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 $$ (3)

donde \( {\hat{n}}^i \) denota vector normal unitario hacia afuera en la interfaz de subdominio Γⁱ, y αⁱ es un parámetro complejo que a menudo puede elegirse como jk₀. Todos los subdominios son luego discretizados por elementos tetraédricos. En cada elemento, expandimos E con funciones base vectoriales N y coeficiente de campo eléctrico desconocido E como

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

donde s denota el número de funciones base vectoriales en cada elemento tetraédrico. s se elige como 6 para función base tradicional de bajo orden basada en borde, mientras que es mayor para función base vectorial de alto orden, ya que se introducen funciones base adicionales basadas en cara o volumen.

Aplicando el método de Galerkin, la ecuación de matriz FEM en Ωⁱ sobre el coeficiente de campo eléctrico desconocido Eⁱ puede obtenerse como

$$ \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) $$ (5)

donde las subíndices notaciones c y r distinguen los grados de libertad (DOFs) de esquina y los DOFs restantes, que extrae los DOFs de esquina como variable primal para construir el esquema dual-primal (DP). Aquí, K es la matriz de sistema FEM y f es el vector de excitación. B es una matriz booleana que extrae los DOFs de interfaz de un subdominio. λ es una variable dual generada de expandir Λ_i, que también se llama multiplicador de Lagrange.

Entonces, los subdominios adyacentes pueden interconectarse forzando la continuidad del campo eléctrico y magnético tangencial en las interfaces. Ensamblamos todas las interfaces de subdominio y eliminamos todos los desconocidos internos de subdominio Eⁱ. Una ecuación de interfaz global reducida sobre la variable dual λ puede obtenerse como

$$ \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 $$ (6)

La Ecuación (6) puede resolverse por métodos iterativos, como el método de residuo mínimo generalizado (GMRES). \( {\tilde{\mathbf{K}}}_{cc} \) es el sistema global de esquina, que puede acelerar la convergencia iterativa en espacio primal. Precondicionador adecuado puede emplearse para mejorar la tasa de convergencia iterativa, como inversa aproximada y descomposición LU incompleta. Una vez que la variable dual λ se resuelve, el campo eléctrico dentro de cada subdominio puede evaluarse independientemente por (5). Para la construcción de la matriz global \( {\tilde{\mathbf{K}}}_{rr} \), se necesita construir la función de Green numérica de subdominio \( {\mathbf{Z}}_{rr}^i \) con una forma de

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

donde la inversa de la matriz FEM de subdominio \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) está incluida. Además, para matrices \( {\tilde{\mathbf{K}}}_{rc} \), \( {\tilde{\mathbf{K}}}_{cr} \), y \( {\tilde{\mathbf{K}}}_{cc} \) y vectores \( {\tilde{f}}_r \) y \( {\tilde{f}}_c \), la \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) se requiere calcular. Las construcciones de \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) en etapa de pre-procesamiento así como sus productos matriz-vector (MVPs) en etapa de solución iterativa son computacionalmente costosos. Aunque \( {\mathbf{K}}_{rr}^i \) es dispersa, \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) son densas, lo que significa altos costos computacionales. A continuación, se introduce un método de esparsificación de bajo rango para acelerar el cómputo de \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \). Dado que algunas sub-matrices en el sistema de interfaz global pueden representarse en forma de matriz de bajo rango, su cómputo puede realizarse con algoritmo de bajo rango, lo que mejora el rendimiento del FETI-DP. Se puede ver que el FETI-DP proporciona operaciones de subdominio independientes, por lo que puede explotar computación paralela para mejorar eficiencia. Para un esquema paralelo eficiente, un principio de división de dominio es hacer el número de DOFs en cada subdominio lo más balanceado posible. Por tanto, el tamaño de subdominios debe relacionarse con la densidad de discretización de malla. Usualmente, subdominios pequeños se adoptan en áreas finamente malladas, mientras que subdominios grandes se adoptan en áreas groseramente malladas.

Esparsificación de bajo rango

Se propone un enfoque de esparsificación de bajo rango para proporcionar una forma de datos dispersos para mejorar la eficiencia FETI-DP. Aquí, datos dispersos significa que estas matrices no son realmente dispersas pero lo son en el sentido de que ciertos sub-bloques de ellas pueden representarse por formas de matriz de descomposición de bajo rango como

$$ \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) $$ (8)

donde X y Y están en formas de matriz completa, y el rango k es mucho menor que m y n. La matriz \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) puede representarse por formas de matriz de datos dispersos ya que posee la propiedad de matriz de un operador integral. Por tanto, siempre que \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) posea propiedad de bajo rango en subdominio dado, puede calcularse y almacenarse eficientemente en formas de datos dispersos con el enfoque de esparsificación de bajo rango, lo que acelera los MVPs en la solución iterativa.

Los procesos del enfoque de esparsificación de bajo rango pueden dividirse en los siguientes pasos: (1) construir un árbol de cluster subdividiendo el conjunto de funciones base en cada subdominio, (2) construir un árbol de cluster de bloques por interacción de dos árboles de cluster, (3) generar una forma de datos dispersos de \( {\mathbf{K}}_{rr}^i \) por una condición de admisibilidad, (4) realizar algoritmos formateados de bajo rango para obtener la representación de datos dispersos de \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}}^{-1} \), y (5) entrar en la solución de sistemas FETI-DP por algoritmo de datos dispersos. Precondicionador adecuado puede emplearse para acelerar la solución. Debe notarse que la factorización LU de datos dispersos \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}}={\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}}{\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}} \) se adopta para reemplazar la inversión de matriz \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}}^{-1} \). Una técnica de disección anidada se emplea paraさらに mejorar la eficiencia de la esparsificación de bajo rango. La disección anidada usa separadores para producir sub-bloques cero fuera de diagonal grandes, que mantendrán cero durante la factorización LU para que pueda reducir significativamente los rellenos.

Para generar \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}} \), primero construimos un árbol de cluster T_I por subdivisión recursiva del conjunto de funciones base basadas en borde de subdominio I = {1,2,……N} usando caja delimitadora. Con la disección anidada, un cluster t dentro de la caja delimitadora correspondiente se divide en tres sucesores {s₁, s_sep, s₂}, donde s₁ y s₂ son los conjuntos de índices de las dos cajas delimitadoras desconectadas y s_sep es el conjunto de índices del separador. La Figura 2a muestra un ejemplo simple de este proceso. Luego, un árbol de cluster de bloques T_{I × I} puede construirse interactuando dos árboles de cluster T_I, como se muestra en Fig. 2b, que puede elegirse como el árbol de cluster del conjunto de funciones base basadas en borde original y el del conjunto de funciones base de prueba en el método de Galerkin. A continuación, necesitamos introducir una condición de admisibilidad basada en la disección anidada para distinguir bloques completos, bloques de descomposición de bajo rango y bloques cero fuera de diagonal en T_{I × I} [23]. Así, \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}} \) puede producirse llenando los bloques correspondientes con las entradas no cero de \( {\mathbf{K}}_{rr}^i \). Finalmente, la factorización LU de datos dispersos de \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}}={\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}}{\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}} \) puede calcularse recursivamente desde

$$ {\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] $$ (9)

Construcciones de un árbol de cluster y un árbol de cluster de bloques de 4 niveles basados en disección anidada. a Construcción de un árbol de cluster por subdivisión recursiva del conjunto de funciones base basadas en borde I = {1,2,…18}. b Construcción de un árbol de cluster de bloques donde bloques blancos son matrices cero y bloques verdes pueden ser matrices completas o matrices de descomposición de bajo rango

donde aritméticas de matriz completa convencionales son reemplazadas por sus contrapartes de datos dispersos [28]. Un error de truncamiento adaptativo ε_t se emplea para controlar la precisión de aproximaciones de bajo rango. Los factores LU obtenidos \( {\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}} \) y \( {\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}} \) se almacenan y usan para construir \( {\mathbf{Z}}_{rr}^i \) por

$$ {\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 $$ (10)

donde \( {\mathbf{B}}_r^i{\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}}^{-1} \) y \( {\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}}^{-1}{\mathbf{B}}_r^i \) pueden calcularse por solucionadores triangulares superior e inferior de datos dispersos. Los \( {\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}} \), \( {\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}} \), y \( {\mathbf{Z}}_{rr}^i \) entran en el cálculo FETI-DP con sustituciones hacia adelante y hacia atrás (FBSs) de datos dispersos y MVPs de datos dispersos.

Fuerza y potencial óptico

Según la teoría de electrodinámica, la fuerza óptica puede evaluarse por el tensor de estrés de Maxwell (MST) que revela la relación entre campo electromagnético y momento mecánico [29]. Una vez obtenida la distribución de campo electromagnético alrededor del objeto, la fuerza óptica puede calcularse integrando MST sobre una superficie cerrada que rodea el objeto. Basado en la distribución de campo eléctrico obtenida, el MST en cualquier coordenada puede construirse por

$$ \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] $$ (11)

donde el superíndice asterisco denota el conjugado de campo eléctrico o magnético, ε son μ son la permitividad y permeabilidad, y \( \overleftrightarrow{\mathbf{I}} \) es una matriz identidad 3 × 3. Por el producto externo de vectores, la forma tensorial de \( \overleftrightarrow{\mathbf{T}} \) puede escribirse como

$$ \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] $$ (12)

donde el subíndice x, y, z denota los componentes en tres direcciones. Según la expansión de E descrita en (4), las entradas de MST T_mn (m, n = x, y, z) pueden convertirse en formas expandidas en el cálculo FETI-DP como

$$ {\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}} $$ (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. $$ (14)

donde ω es la frecuencia angular; N y s han sido descritos en Ec. (4).

Finalmente, la fuerza óptica ejercida sobre el objeto puede calcularse integrando el MST sobre cualquier superficie cerrada que lo rodee por

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

Notar que el cálculo de fuerza óptica también puede implementarse en paralelo, ya que la integral del MST se asigna a subdominios correspondientes. Para un atrapamiento óptico estable, una de las condiciones principales es que la fuerza de gradiente debe ser mayor que la fuerza de dispersión. En otras palabras, la dirección de la fuerza total debe ser idéntica a la de la fuerza de gradiente, que siempre apunta a la posición donde la intensidad de campo eléctrico es más fuerte.

El potencial óptico es otro parámetro atractivo revelando la estabilidad del atrapamiento óptico. Basado en la fuerza óptica obtenida, la profundidad de potencial óptico U en posición r₀ puede calcularse por

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

donde el subíndice ∞ denota infinito definido como punto de referencia con potencial cero. El valor de U puede representarse por k_BT, donde k_B denota la constante de Boltzmann de 1,3806488 × 10⁻²³ J/K y T es la temperatura ambiente. Generalmente, la partícula puede superar el movimiento browniano en solución y ser atrapada establemente cuando U > 1 k_BT se satisface. De lo contrario, la partícula no puede ser atrapada establemente. Dado que la fuerza óptica total incluye el componente conservativo de fuerza de gradiente y el componente no conservativo de fuerza de dispersión, la fuerza óptica total F de (15) es no conservativa [30, 31]. Sin embargo, siempre que el movimiento de la nanopartícula esté restringido a una dimensión, esto produce una definición inequívoca de un potencial óptico de (16), aunque la fuerza óptica total sea no conservativa.

Resultados y discusión

Se presentan tres ejemplos para demostrar la efectividad del método propuesto. Dado que metales nobles se usan comúnmente para excitar el plasmon de superficie, seleccionamos materiales representativos de oro y plata para los análisis. El primer ejemplo calcula la fuerza óptica de nanopartícula de plata para verificar la precisión del método FETI-DP propuesto. El segundo y tercer ejemplo simulan y discuten el atrapamiento óptico de nanopartículas de oro. Para todos los ejemplos, el dominio infinito se trunca con ABC, y las distancias entre el ABC y los objetos se establecen en una longitud de onda, lo cual es suficiente para lograr una precisión aceptable. Todos los cálculos se realizan en una estación de trabajo Dell equipada con procesadores Intel Xeon de 3,6 GHz.

Nanocápsula de plata

Primero se considera un objeto de nanocápsula de plata para probar la precisión y eficiencia del método FETI-DP propuesto en la predicción de fuerza óptica. La Figura 3a y b presenta su configuración y dimensiones. Los parámetros constitutivos de plata son todos valores medidos tomados de [32]. Para implementar el esquema FETI-DP, el dominio de análisis completo se divide primero en 24 subdominios. Se requieren mallas más densas cerca de la superficie metálica para modelar el efecto de mejora de campo local plasmonico. Se adoptan elementos tetraédricos para la discretización, lo que lleva a totalmente 6,9 × 10⁵ incógnitas, incluyendo 4,1 × 10⁴ incógnitas duales y 313 incógnitas de esquina. La luz incide a lo largo de la dirección +z, mientras que la dirección de polarización del campo eléctrico es -x.

Configuración de una estructura de nanocápsula de plata. a Vista 3D. b Vista frontal y dimensiones, donde R = 30 nm y h = 60 nm

Primero, cambiamos la longitud de onda de luz incidente λ de 200 nm a 400 nm para simular las fuerzas ópticas ejercidas sobre la nanocápsula. Dado que el FETI-DP trabaja en dominio de frecuencia, las fuerzas ópticas se calculan en 15 puntos de muestreo de frecuencia. La Figura 4 muestra la curva calculada de fuerzas ópticas ejercidas sobre la nanocápsula de plata. Para indicar la precisión del FETI-DP, los resultados de fuerza óptica del FETI-DP se comparan con los del software comercial Lumerical FDTD Solutions [33], y se puede observar buen acuerdo.

Resultados de las fuerzas ópticas ejercidas sobre la nanocápsula de plata, variando con la longitud de onda λ de luz incidente, incluyendo los resultados del FETI-DP y las soluciones FDTD del software comercial

Luego, se prueba el rendimiento del FETI-DP para diferentes números de subdominios. Aumentamos el número de subdominios de 4 a 24 manteniendo la densidad de discretización. Asignamos cada procesador para lidiar con un subdominio. La Tabla 1 reporta el tiempo usado para la construcción de la ecuación de interfaz global Ec. (6) y el tiempo total de solución. Se puede ver que el FETI-DP puede explotar completamente recursos de computación paralela y mejorar significativamente la eficiencia de solución. Además, la precisión del FETI-DP con el número de subdominios aumentando también se examina y reporta en Tabla 1. Aquí, la precisión se define por el error relativo de 2-norma de la fuerza óptica como δ_OF = ‖OF_i − OF_ref‖/‖OF_ref‖, donde OF_i es la fuerza óptica usando i subdominios y OF_ref denota la fuerza óptica de referencia usando dos subdominios. Se puede ver que la precisión se mantiene casi constante con el número de subdominios aumentando.

Dímero de nanoesfera de oro

El segundo ejemplo analiza el atrapamiento óptico de una nanoesfera de oro usando un dímero de nanoesfera de oro. Los efectos plasmonicos en el hueco del dímero pueden mejorar efectivamente la fuerza óptica para atrapar nanopartícula. La Figura 5a y b da la configuración y dimensiones de este sistema. Los parámetros constitutivos de oro son todos valores medidos tomados de [32]. El medio circundante es agua con un índice de refracción relativo de n = 1,33. La luz incidente es una onda plana con potencia de 10 mW/µm², la dirección de polarización del campo eléctrico es +x, y la dirección de incidencia es -z. La fuerza óptica ejercida sobre el objeto nanoesfera se calcula por el método FETI-DP. Para la implementación FETI-DP, el dominio computacional completo se divide en 32 subdominios y discretizado por mallas tetraédricas, lo que resulta en 3,5 × 10⁶ incógnitas, incluyendo 1,6 × 10⁵ incógnitas duales y 1738 incógnitas de esquina.

Configuración de un sistema de atrapamiento óptico de un dímero de nanoesfera de oro en agua. a Vista 3D. b Vista frontal y dimensiones, donde R = 25 nm, r = 5 nm, y g = 2 nm

Primero, probamos el rendimiento paralelo del FETI-DP propuesto usando varios números de procesadores. La Tabla 2 reporta el tiempo de solución para Ec. (6) así como el tiempo total de solución. Además, los speedups para la computación paralela también se proporcionan en Tabla 2. Aquí, el speedup se define por

$$ \mathrm{Speed}\ \mathrm{up}=\raisebox{1ex}{${T}_1$}\!\left/ \!\raisebox{-1ex}{${T}_{N_p}$}\right. $$ (17)

donde \( {T}_{N_p} \) denota el tiempo total de reloj de pared usando N_p procesadores. Se puede ver que el FETI-DP mejora significativamente la eficiencia de solución y exhibe buen speedup paralelo. Para este gran número de incógnitas, el uso total de memoria de todos los procesadores es solo 57,2 GB.

Luego, se examina la efectividad del enfoque de esparsificación de bajo rango. Con la esparsificación de bajo rango, la matriz de subdominio puede factorizarse por algoritmo de datos dispersos y almacenarse como matrices de datos dispersos. El tiempo de construcción y uso de memoria son solo 18 s y 0,5 GB, mientras que son 67 s y 1,7 GB por algoritmo de matriz convencional. Se puede ver que obtenemos 72 % de ahorro de tiempo y 70 % de compresión de memoria. Relacionado con el uso de memoria, los MVPs subsecuentes también pueden obtener 70 % de ahorro de tiempo.

Luego, el FETI-DP se prueba para el cálculo de fuerza óptica con la longitud de onda λ variando de 277 nm a 818 nm. En la práctica, los análisis de fuerza óptica bajo luz incidente de diferentes longitudes de onda son a menudo necesarios para buscar la longitud de onda de resonancia plasmonica, donde ocurre mejora drástica de campo y se puede obtener la fuerza óptica más fuerte. Se consideran dos casos con la nanoesfera ubicada en (0, 0, 20 nm) y (0, 0, -20 nm). La Figura 6a y b traza las fuerzas ópticas calculadas ejercidas sobre la nanoesfera para diferentes λ. Se puede ver que la fuerza óptica máxima ocurre en λ = 472 nm, que es la longitud de onda de resonancia plasmonica. La fuerza óptica en esta longitud de onda de resonancia se mejora casi 40 veces contra la de longitud de onda no resonante. Además, la fuerza óptica siempre apunta al hueco del dímero, como se muestra en Fig. 6, donde la intensidad de campo eléctrico es más fuerte. También es la dirección de fuerza de gradiente para atrapar el objeto. La Figura 7a y b muestra las distribuciones calculadas de mejora de campo eléctrico en la longitud de onda no resonante de λ = 300 nm y la longitud de onda de resonancia de λ = 472 nm, respectivamente. Se puede ver que la intensidad de campo eléctrico ha aumentado casi 250 veces debido al efecto de resonancia plasmonica.

Resultados calculados de fuerzas ópticas ejercidas sobre la nanoesfera en el sistema de dímero de nanoesfera de oro, variando con la longitud de onda λ de luz incidente. a La nanoesfera objeto está ubicada en (0, 0, 20 nm). b La nanoesfera objeto está ubicada en (0, 0, -20 nm)

Distribuciones de mejora de campo eléctrico en el plano xoz para el sistema de dímero de nanoesfera de oro. a λ = 300 nm (longitud de onda no resonante). b λ = 472 nm (longitud de onda de resonancia)

Además, la fuerza óptica y potencial óptico se calculan con la nanoesfera moviéndose de (0, 0, -30 nm) a (0, 0, -17 nm) a lo largo del eje z. Dado que el comportamiento más típico e interesante de fuerzas y potenciales de atrapamiento son aquellos actuando a lo largo de la dirección z, aquí consideramos el potencial de atrapamiento axial por integración a lo largo del eje z. Porque el movimiento de la nanopartícula está restringido a una dimensión, la definición de un potencial óptico es inequívoca de (16), aunque la fuerza óptica total de (15) sea no conservativa. Como se muestra en Fig. 8a, b, con la nanoesfera moviéndose al hueco del dímero, la fuerza óptica y profundidad de potencial óptico aumentan obviamente. En la posición de (0, 0, -17 nm), se produce una profundidad de potencial óptico de 4,6 k_BT, que es suficiente para superar el movimiento browniano en agua para lograr atrapamiento óptico estable.

Las fuerzas ópticas y potenciales ópticos ejercidos sobre la nanoesfera en el sistema de dímero de nanoesfera de oro, cuando la nanoesfera se mueve de (0, 0, -30 nm) a (0, 0, -17 nm). a Las fuerzas ópticas. b Los potenciales ópticos

Finalmente, probamos los efectos del sustrato dieléctrico para este ejemplo. Las fuerzas ópticas se calculan con y sin sustrato, respectivamente. Para ambos casos, la nanoesfera está ubicada en (0, 0, -20 nm) y la longitud de onda de incidencia se elige como la longitud de onda de resonancia. Para el caso sin sustrato, el resultado calculado de la fuerza óptica es |F₀| = 0,769 pN. Para el caso con sustrato, el dímero de nanoesfera de oro se coloca en un sustrato dieléctrico con espesor de 60 nm y permitividad relativa de ε_r = 2,25. El resultado calculado de la fuerza óptica es |F₁| = 0,761 pN. El error relativo entre estos dos resultados de fuerzas ópticas es约 1,0 × 10⁻², definido como |F₁ − F₀|/|F₀|.

Dímero de cono truncado de oro

El tercer ejemplo trata el atrapamiento óptico de una nanoesfera de oro usando un dímero de cono truncado de oro. La Figura 9 da la configuración y dimensiones de este sistema. Los parámetros constitutivos de oro se toman de [32]. El sustrato dieléctrico tiene una permitividad relativa de ε_r = 2,25. El medio circundante es agua con un índice de refracción relativo de n = 1,33. La luz incidente es onda plana con potencia de 10 mW/µm², la dirección de polarización del campo eléctrico es +x, y la dirección de incidencia es -z. El dominio computacional completo se divide en 32 subdominios y discretizado por mallas tetraédricas, lo que lleva a 3,1 × 10⁶ incógnitas, incluyendo 1,3 × 10⁵ incógnitas duales y 1227 incógnitas de esquina.

Configuración de un sistema de atrapamiento óptico de un dímero de cono truncado de oro basado en un sustrato dieléctrico en agua. a Vista 3D. b Vista frontal y dimensiones, donde UR = 20 nm, LR = 30 nm, h = 35 nm, y g = 2 nm

Primero, analizamos las fuerzas ópticas cambiando λ de 277 nm a 818 nm. La Figura 10 traza las fuerzas ópticas calculadas ejercidas sobre la nanoesfera para diferentes λ. La nanoesfera está ubicada en (0, 0, 35 nm). Se puede ver que la fuerza óptica máxima ocurre en λ = 464 nm, que es la longitud de onda de resonancia plasmonica, y la fuerza óptica aquí se mejora casi 30 veces en longitud de onda no resonante. Además, la fuerza óptica total siempre apunta a -z, como se muestra en Fig. 10, que es la dirección de la fuerza de gradiente. Esto confirma que la fuerza de gradiente es mayor que la fuerza de dispersión, que es una de las condiciones para que la nanoesfera pueda ser atrapada establemente. La Figura 11a y b presenta las distribuciones calculadas de campo eléctrico en la longitud de onda no resonante de λ = 300 nm y la longitud de onda de resonancia de λ = 464 nm, respectivamente. Se puede ver que la intensidad de campo eléctrico ha aumentado casi 500 veces debido a la resonancia de plasmon de superficie localizada.

Resultados calculados de fuerzas ópticas ejercidas sobre la nanoesfera en el sistema de dímero de cono truncado de oro, variando con λ. La nanoesfera está ubicada en (0, 0, 35 nm)

Distribuciones de mejora de campo eléctrico en el plano xoz para el sistema de dímero de cono truncado de oro. a λ = 300 nm (longitud de onda no resonante). b λ = 464 nm (longitud de onda de resonancia)

Luego, la ubicación de la nanoesfera se cambia a 0, 5 y 35 nm para observar la fuerza óptica. La Figura 12 da las fuerzas ópticas calculadas ejercidas sobre la nanoesfera, donde se puede observar un componente y obvio de fuerza óptica, mientras existe un componente z mayor de fuerza óptica. La fuerza óptica total aún apunta a la posición con el campo eléctrico más fuerte para atrapar la nanoesfera.

Resultados calculados de fuerzas ópticas ejercidas sobre la nanoesfera en el sistema de dímero de cono truncado de oro variando λ. La nanoesfera está ubicada en (0, 5 nm, 35 nm)

Además, analizamos el potencial óptico con la nanoesfera moviéndose de (0, 0, 50 nm) a (0, 0, 20 nm) a lo largo del eje z. Aquí, consideramos el potencial de atrapamiento axial a lo largo de la dirección z, que restringe el movimiento de la nanopartícula a una dimensión y lleva a una definición inequívoca de potencial óptico. Tanto la fuerza óptica como el potencial se calculan. Como se puede observar en Fig. 13a, b, con la nanoesfera moviéndose al hueco del dímero, la fuerza óptica y la profundidad de potencial óptico aumentan obviamente. En (0, 0, 20 nm), se obtiene una profundidad de potencial óptico de 3,8 k_BT, que es suficiente para superar el movimiento browniano en agua para lograr atrapamiento óptico estable.

Las fuerzas ópticas y potenciales ópticos ejercidos sobre la nanoesfera en el sistema de dímero de cono truncado de oro, cuando la nanoesfera se mueve de (0, 0, 50 nm) a (0, 0, 20 nm). a Las fuerzas ópticas. b Los potenciales ópticos

Finalmente, probamos los costos computacionales del FETI-DP cambiando el número de incógnitas de 1,0 millón a 3,2 millones basado en diferentes tamaños de malla. En la práctica, las pruebas bajo diferentes densidades de malla son usualmente necesarias para cumplir diferentes requisitos de precisión. Un problema complejo de tan gran escala trae grandes desafíos a métodos numéricos convencionales. Sin embargo, el FETI-DP puede manejar fácilmente este problema. Treinta y dos procesadores se emplean para la simulación FETI-DP, mientras cada procesador lidia con un subdominio. La Tabla 3 reporta los costos computacionales del FETI-DP. Se puede ver que el FETI-DP exhibe alta eficiencia de simulación y bajo requerimiento de memoria.

Conclusión

Se propone un método FETI-DP combinado con esparsificación de bajo rango para la predicción y análisis de atrapamiento óptico de nanopartículas metálicas. El método propuesto proporciona problemas de subdominio completamente desacoplados, que convierte un problema complejo de gran escala en una serie de problemas simples de pequeña escala. Es adecuado para computación paralela y puede mejorar significativamente la eficiencia de simulación numérica. Los ejemplos demuestran que el método propuesto exhibe excelente rendimiento de computación de gran escala y es adecuado para la simulación rápida y precisa de atrapamiento óptico a nanoescala.

Disponibilidad de datos y materiales

Todos los datos generados o analizados durante este estudio están incluidos en este artículo.

Abreviaturas

ABC:

Condición de frontera absorbente

DOF:

Grados de libertad

FDTD:

Diferencias finitas en dominio temporal

FEM:

Método de elementos finitos

FETI-DP:

Desgarro e interconexión dual-primal de elementos finitos

MST:

Tensor de estrés de Maxwell

MVP:

Producto matriz-vector

Correcciones a las fórmulas de números mágicos para nanoclústeres cúbicos centrados en el cuerpo y nuevos datos de cubo truncado Nanohelice de doble compuerta como superred binaria ajustable: Ingeniería de bandas y perspectivas optoelectrónicas