Simulation précise et efficace du piégeage optique nanoscopique par FETI-DP et tenseur de Maxwell

Résumé

La simulation numérique joue un rôle important pour la prédiction du piégeage optique basé sur des pinces optiques nano-plasmoniques. Cependant, les structures complexes et l'amélioration drastique du champ local des effets plasmoniques posent de grands défis aux méthodes numériques traditionnelles. Dans cet article, une méthode de simulation numérique précise et efficace basée sur une méthode de déchirement et d'interconnexion d'éléments finis dual-primal (FETI-DP) et le tenseur de contrainte de Maxwell est proposée, pour calculer la force et le potentiel optiques pour le piégeage de nanoparticules. Une approche de parcimonie de bas rang est introduite pour améliorer davantage les performances de simulation FETI-DP. La méthode proposée peut décomposer un problème complexe à grande échelle en problèmes simples à petite échelle en utilisant une division de domaine non chevauchante et une discrétisation de maillage flexible, ce qui présente une haute efficacité et parallélisabilité. Les résultats numériques montrent l'efficacité de la méthode proposée pour la prédiction et l'analyse du piégeage optique à l'échelle nanométrique.

Introduction

Les pinces optiques plasmoniques basées sur les plasmons de surface (PS) attirent beaucoup d'attention et ont été largement appliquées pour capturer des nanoparticules [1-6]. Le PS est un phénomène de résonance causé par le couplage de la lumière incidente d'une longueur d'onde spécifique et des électrons libres à l'interface des métaux et diélectriques [7]. Les PS permettent aux pinces optiques de franchir la limite de diffraction. De plus, l'amélioration drastique du champ local des PS peut réduire la demande d'intensité de la lumière incidente [7, 8]. Cependant, les PS sont étroitement liés au matériau et aux dimensions des objets ainsi qu'à la longueur d'onde de la lumière incidente, ce qui nécessite un grand nombre d'expériences pour déterminer les paramètres optimaux des pinces optiques PS en pratique. Basé sur cela, la méthode de simulation joue un rôle de plus en plus important comme moyen auxiliaire pour la conception et l'optimisation des pinces optiques PS [9]. Dans ces simulations, le calcul de la force optique est requis pour prédire un piégeage stable. Pour des objets réguliers tels que les sphères, la force optique peut être dérivée analytiquement de la théorie de Lorenz-Mie généralisée [10, 11]. Cependant, pour des objets avec des configurations complexes, des méthodes numériques qui résolvent rigoureusement les équations de Maxwell gouvernantes sont nécessaires pour modéliser les champs électromagnétiques et la force et le potentiel optiques subséquents.

Ces méthodes numériques peuvent être principalement catégorisées en méthodes d'équations différentielles (MED) et méthodes d'équations intégrales (MEI) [12-15]. Comparées aux MEI, les méthodes d'équations différentielles (MED) montrent des capacités supérieures pour traiter les géométries et composants complexes. Les MED ont aussi l'avantage d'un calcul direct de la distribution de champ proche, qui joue un rôle important dans l'analyse des PS. En tant que MED représentative, la méthode de différences finies dans le domaine temporel (FDTD) est implémentée dans le domaine temporel, ce qui peut facilement obtenir des informations à large bande et des réponses transitoires [16, 17]. Cependant, le FDTD exige un modèle dispersif précis pour décrire les propriétés matérielles dépendantes de la fréquence dans les PS, tandis que la précision de la solution FDTD dépend fortement de la précision d'approximation de ce modèle dispersif [18]. De plus, le FDTD repose sur des maillages structurés, qui mènent souvent à une erreur d'escalier pour les surfaces courbes. En tant qu'autre MED représentative, la méthode des éléments finis (MEF) a été largement adoptée car elle peut facilement gérer les matériaux dispersifs dans le domaine fréquentiel et éliminer l'erreur d'escalier par maillage non structuré [19-22]. Comparé au FDTD, le MEF peut directement adopter les paramètres matériels mesurés sans aucun modèle dispersif approximatif. Cependant, les améliorations drastiques du champ local dans les PS nécessitent des maillages fins dans la discrétisation MEF. De plus, les objets avec de grandes dimensions et de multiples objets augmenteront dramatiquement le nombre d'inconnues. Ces facteurs causeront des systèmes de matrices mal conditionnés et d'énormes consommations computationnelles, ce qui pose de grands défis au MEF traditionnel pour l'analyse du piégeage optique amélioré par PS.

Dans cet article, une méthode FETI-DP (dual-primal finite element tearing and interconnecting) efficace est introduite pour simuler le piégeage optique à l'échelle nanométrique. Le FETI-DP adopte un schéma de décomposition de domaine non chevauchant, qui divise un problème complexe original à grande échelle en une série de problèmes simples à petite échelle pour les conquérir. Il impose une condition de transmission aux interfaces de sous-domaines pour assurer la continuité du champ, et introduit une variable duale pour réduire le problème original tridimensionnel (3D) en un problème bidimensionnel (2D) par multiplicateur de Lagrange. Les variables primales aux coins des sous-domaines sont extraites pour accélérer le taux de convergence de la solution itérative du problème dual [23-26]. Une approche de parcimonie de bas rang est développée pour améliorer les performances du FETI-DP. Elle utilise des algorithmes parcimonieux en données pour améliorer l'efficacité pour résoudre les problèmes de sous-domaines et le problème dual [27, 28]. La méthode proposée fournit des sous-domaines pleinement découplés, ce qui permet la simulation parallèle de la force optique pour le piégeage de nanoparticules. Le tenseur de contrainte de Maxwell (TCM) qui révèle la relation entre le champ électromagnétique et le momentum mécanique est adopté pour évaluer la force optique [29]. Basé sur la force optique obtenue, le potentiel optique peut être calculé pour l'analyse d'un piégeage stable. Comparé aux MEI, la méthode proposée est plus puissante pour traiter les matériaux composés et résoudre le champ proche pour le piégeage optique basé sur PS. Comparé au FDTD, la méthode proposée peut gérer avec précision le matériau métallique dispersif dans les systèmes de piégeage optique basés sur PS et éliminer l'erreur d'escalier pour les objets à frontière courbe. Comparé au MEF, la méthode proposée est appropriée pour le calcul à grande échelle du piégeage optique. Plusieurs exemples sont analysés et les résultats numériques démontrent la précision et l'efficacité de la méthode proposée pour la prédiction et l'analyse du piégeage optique à l'échelle nanométrique.

Méthodes

Formulations FETI-DP

Pour l'implémentation FETI-DP, le domaine computationnel original Ω est d'abord déchiré en une série de sous-domaines non chevauchants Ωⁱ (i = 1, 2, 3⋯, N_s), comme montré en Fig. 1. Dans chaque sous-domaine Ωⁱ, un système d'éléments finis de sous-domaine peut être dérivé de l'équation d'onde vectorielle comme :

$$ \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{dans}\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{sur}\kern0.24em {\Gamma}_{\mathrm{ABC}}^i $$ (2)

Schéma de division de domaine avec sous-domaines non chevauchants dans la méthode FETI-DP. a Domaine original. b Sous-domaines divisés et maillages discrétisés.

où Eⁱ désigne le champ électrique inconnu à résoudre dans \( {\Omega}^i \), k₀ et η₀ sont le nombre d'onde et l'impédance intrinsèque du vide, respectivement, et \( {\mathbf{J}}_i^{\mathrm{imp}} \) est le courant imprimé. \( {\Gamma}_{\mathrm{ABC}}^i \) signifie la condition de frontière absorbante (ABC) pour tronquer la région ouverte infinie. Il doit être noté que k₀ devrait être remplacé par l'impédance d'onde dans le milieu si le milieu environnant n'est pas le vide, ce qui est commun pour le piégeage optique. À l'interface de sous-domaine Γⁱ, une condition de frontière supposée est requise pour générer un problème de valeur de frontière complet dans Ωⁱ. Ici, une condition de transmission de type Robin avec une variable auxiliaire inconnue Λ_i est imposée comme :

$$ {\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{sur}\kern0.36em {\Gamma}^i $$ (3)

où \( {\hat{n}}^i \) désigne le vecteur normal unitaire sortant à l'interface de sous-domaine Γⁱ, et αⁱ est un paramètre complexe qui peut souvent être choisi comme jk₀. Tous les sous-domaines sont ensuite discrétisés par éléments tétraédriques. Dans chaque élément, nous développons E avec des fonctions de base vectorielles N et un coefficient de champ électrique inconnu E comme :

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

où s désigne le nombre de fonctions de base vectorielles dans chaque élément tétraédrique. s est choisi à 6 pour la fonction de base vectorielle d'ordre bas traditionnel basée sur l'arête, tandis qu'il est plus grand pour la fonction de base vectorielle d'ordre élevé, puisque des fonctions de base additionnelles basées sur la face ou le volume sont introduites.

En appliquant la méthode de Galerkin, l'équation matricielle MEF dans Ωⁱ sur le coefficient de champ électrique inconnu Eⁱ peut être obtenue comme :

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

où les notations d'indice c et r distinguent les degrés de liberté (DDL) de coin et les DDL restants, ce qui extrait les DDL de coin comme variable primale pour construire le schéma dual-primal (DP). Ici, K est la matrice de système MEF et f est le vecteur d'excitation. B est une matrice booléenne qui extrait les DDL d'interface d'un sous-domaine. λ est une variable duale générée en développant Λ_i, aussi appelée multiplicateur de Lagrange.

Ensuite, les sous-domaines adjacents peuvent être interconnectés en imposant la continuité du champ électrique et magnétique tangentiel aux interfaces. Nous assemblons toutes les interfaces de sous-domaines et éliminons tous les inconnus internes de sous-domaine Eⁱ. Une équation d'interface globale réduite sur la variable duale λ peut être obtenue comme :

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

L'équation (6) peut être résolue par des méthodes itératives, telles que la méthode des résidus minimaux généralisés (GMRES). \( {\tilde{\mathbf{K}}}_{cc} \) est le système global de coins, qui peut accélérer la convergence itérative dans l'espace primal. Un préconditionneur approprié peut être employé pour améliorer le taux de convergence itératif, tel que l'inverse approché et la décomposition LU incomplète. Une fois la variable duale λ résolue, le champ électrique à l'intérieur de chaque sous-domaine peut être évalué indépendamment par (5). Pour la construction de la matrice globale \( {\tilde{\mathbf{K}}}_{rr} \), il faut construire la fonction de Green numérique de sous-domaine \( {\mathbf{Z}}_{rr}^i \) sous la forme :

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

où l'inverse de la matrice MEF de sous-domaine \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) est inclus. De plus, pour les matrices \( {\tilde{\mathbf{K}}}_{rc} \), \( {\tilde{\mathbf{K}}}_{cr} \), et \( {\tilde{\mathbf{K}}}_{cc} \) et vecteurs \( {\tilde{f}}_r \) et \( {\tilde{f}}_c \), le \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) est requis pour être calculé. Les constructions de \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) à l'étape de pré-traitement ainsi que leurs produits matrice-vecteur (PMV) à l'étape de solution itérative sont coûteux computationnellement. Bien que \( {\mathbf{K}}_{rr}^i \) soit parcimonieuse, \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) sont denses, ce qui signifie des coûts computationnels élevés. Ensuite, une méthode de parcimonie de bas rang est introduite pour accélérer le calcul de \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \). Puisque certaines sous-matrices dans le système d'interface global peuvent être représentées sous forme de matrice de bas rang, leur calcul peut être effectué avec un algorithme de bas rang, ce qui améliore les performances du FETI-DP. On peut voir que le FETI-DP fournit des opérations de sous-domaine indépendantes, de sorte qu'il peut exploiter le calcul parallèle pour améliorer l'efficacité. Pour un schéma parallèle efficace, un principe de division de domaine est de rendre le nombre de DDL dans chaque sous-domaine aussi équilibré que possible. Donc, la taille des sous-domaines devrait se rapporter à la densité de discrétisation de maillage. Habituellement, de petits sous-domaines sont adoptés dans les zones finement maillées, tandis que de grands sous-domaines sont adoptés dans les zones grossièrement maillées.

Parcimonie de bas rang

Une approche de parcimonie de bas rang est proposée pour fournir une voie parcimonieuse en données pour améliorer l'efficacité du FETI-DP. Ici, parcimonieux en données signifie que ces matrices ne sont pas réellement parcimonieuses mais elles sont parcimonieuses dans le sens que certains sous-blocs d'entre elles peuvent être représentés par des formes de matrices de décomposition de bas rang comme :

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

où X et Y sont sous formes matricielles pleines, et le rang k est beaucoup plus petit que m et n. La matrice \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) peut être représentée par des formes matricielles parcimonieuses en données car elle possède la propriété matricielle d'un opérateur intégral. Donc, pourvu que \( {\left({\mathbf{K}}_{rr}^i\right)}^{-1} \) possède la propriété de bas rang dans un sous-domaine donné, elle peut être calculée et stockée efficacement sous formes parcimonieuses en données avec l'approche de parcimonie de bas rang, ce qui accélère les PMV dans la solution itérative.

Les processus de l'approche de parcimonie de bas rang peuvent être divisés en les étapes suivantes : (1) construire un arbre de clusters en subdivisant l'ensemble de fonctions de base dans chaque sous-domaine, (2) construire un arbre de clusters de blocs par interaction de deux arbres de clusters, (3) générer une forme parcimonieuse en données de \( {\mathbf{K}}_{rr}^i \) par une condition d'admissibilité, (4) effectuer des algorithmes formatés de bas rang pour obtenir la représentation parcimonieuse en données de \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}}^{-1} \), et (5) entrer la solution des systèmes FETI-DP par algorithme parcimonieux en données. Un préconditionneur approprié peut être employé pour accélérer la solution. Il doit être noté que la factorisation LU parcimonieuse en données \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}}={\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}}{\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}} \) est adoptée pour remplacer l'inversion de matrice \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}}^{-1} \). Une technique de dissection emboîtée est employée pour améliorer davantage l'efficacité de la parcimonie de bas rang. La dissection emboîtée utilise des séparateurs pour produire de grands sous-blocs hors diagonale nuls, qui resteront nuls pendant la factorisation LU de sorte qu'elle peut réduire significativement les remplissages.

Pour générer \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}} \), nous construisons d'abord un arbre de clusters T_I par subdivision récursive de l'ensemble de fonctions de base basées sur l'arête du sous-domaine I = {1,2,……N} en utilisant une boîte englobante. Avec la dissection emboîtée, un cluster t dans la boîte englobante correspondante est divisé en trois successeurs {s₁, s_sep, s₂}, où s₁ et s₂ sont les ensembles d'indices des deux boîtes englobantes déconnectées et s_sep est l'ensemble d'indices du séparateur. La Figure 2a montre un exemple simple de ce processus. Ensuite, un arbre de clusters de blocs T_{I × I} peut être construit en interagissant deux arbres de clusters T_I, comme montré en Fig. 2b, qui peut être choisi comme l'arbre de clusters de l'ensemble de fonctions de base basées sur l'arête original et celui de l'ensemble de fonctions de base de test dans la méthode de Galerkin. Ensuite, nous devons introduire une condition d'admissibilité basée sur la dissection emboîtée pour distinguer les blocs pleins, les blocs de décomposition de bas rang et les blocs hors diagonale nuls dans T_{I × I} [23]. Ainsi, \( {\left({\mathbf{K}}_{rr}^i\right)}_{\mathrm{DS}} \) peut être produit en remplissant les blocs correspondants avec les entrées non nulles de \( {\mathbf{K}}_{rr}^i \). Finalement, la factorisation LU parcimonieuse en données 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}} \) peut être calculée récursivement depuis :

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

Constructions d'un arbre de clusters et d'un arbre de clusters de blocs de 4 niveaux basés sur la dissection emboîtée. a Construction d'un arbre de clusters par subdivision récursive de l'ensemble de fonctions de base basées sur l'arête I = {1,2,…18}. b Construction d'un arbre de clusters de blocs où les blocs blancs sont des matrices nulles et les blocs verts peuvent être des matrices pleines ou des matrices de décomposition de bas rang.

où les arithmétiques de matrices pleines conventionnelles sont remplacées par leurs homologues parcimonieux en données [28]. Une erreur de troncature adaptative ε_t est employée pour contrôler la précision des approximations de bas rang. Les facteurs LU obtenus \( {\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}} \) et \( {\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}} \) sont stockés et utilisés pour construire \( {\mathbf{Z}}_{rr}^i \) par :

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

où \( {\mathbf{B}}_r^i{\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}}^{-1} \) et \( {\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}}^{-1}{\mathbf{B}}_r^i \) peuvent être calculés par solveurs triangulaires supérieurs et inférieurs parcimonieux en données. Le \( {\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}} \), \( {\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}} \), et \( {\mathbf{Z}}_{rr}^i \) entrent dans le calcul FETI-DP avec substitutions avant et arrière parcimonieuses en données (FBS) et PMV parcimonieux en données.

Force et potentiel optiques

Selon la théorie de l'électrodynamique, la force optique peut être évaluée par le tenseur de contrainte de Maxwell (TCM) qui révèle la relation entre le champ électromagnétique et le momentum mécanique [29]. Une fois la distribution de champ électromagnétique autour de l'objet obtenue, la force optique peut être calculée en intégrant le TCM sur une surface fermée entourant l'objet. Basé sur la distribution de champ électrique obtenue, le TCM à n'importe quelles coordonnées peut être construit par :

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

où l'exposant astérisque dénote le conjugué du champ électrique ou magnétique, ε et μ sont la permittivité et la perméabilité, et \( \overleftrightarrow{\mathbf{I}} \) est une matrice identité 3 × 3. Par le produit externe de vecteurs, la forme tensorielle de \( \overleftrightarrow{\mathbf{T}} \) peut être écrite comme :

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

où l'indice x, y, z dénote les composantes dans trois directions. Selon le développement de E décrit en (4), les entrées du TCM T_mn (m, n = x, y, z) peuvent être converties en formes développées dans le calcul FETI-DP comme :

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

où ω est la fréquence angulaire ; N et s ont été décrits dans l'Éq. (4).

Finalement, la force optique exercée sur l'objet peut être calculée en intégrant le TCM sur toute surface fermée l'entourant par :

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

Noter que le calcul de la force optique peut aussi être implémenté en parallèle, puisque l'intégrale du TCM est assignée aux sous-domaines correspondants. Pour un piégeage optique stable, une des conditions principales est que la force de gradient devrait être supérieure à la force de diffusion. En d'autres termes, la direction de la force totale devrait être identique à celle de la force de gradient, qui pointe toujours vers la position où l'intensité du champ électrique est la plus forte.

Le potentiel optique est un autre paramètre attractif révélant la stabilité du piégeage optique. Basé sur la force optique obtenue, la profondeur de potentiel optique U à la position r₀ peut être calculée par :

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

où l'indice ∞ dénote l'infini défini comme point de référence avec potentiel nul. La valeur de U peut être représentée par k_BT, où k_B dénote la constante de Boltzmann de 1,3806488 × 10⁻²³ J/K et T est la température ambiante. Généralement, la particule peut surmonter le mouvement brownien en solution et être piégée stablement quand U > 1 k_BT est satisfait. Sinon, la particule ne peut pas être piégée stablement. Puisque la force optique totale inclut la composante de force de gradient conservative et la composante de force de diffusion non-conservative, la force optique totale F de (15) est non-conservative [30, 31]. Cependant, pourvu que le mouvement de la nanoparticule soit restreint à une dimension, cela produit une définition non ambiguë d'un potentiel optique de (16), même si la force optique totale est non-conservative.

Résultats et discussion

Trois exemples sont présentés pour démontrer l'efficacité de la méthode proposée. Puisque les métaux nobles sont couramment utilisés pour exciter le plasmon de surface, nous sélectionnons des matériaux représentatifs or et argent pour les analyses. Le premier exemple calcule la force optique de nanoparticule d'argent pour vérifier la précision de la méthode FETI-DP proposée. Le second et troisième exemples simulent et discutent le piégeage optique de nanoparticules d'or. Pour tous les exemples, le domaine infini est tronqué avec ABC, et les distances entre l'ABC et les objets sont fixées à une longueur d'onde, ce qui est suffisant pour atteindre une précision acceptable. Tous les calculs sont effectués sur une station de travail Dell équipée de processeurs Intel Xeon 3,6 GHz.

Nanocapsule d'argent

Un objet nanocapsule d'argent est d'abord considéré pour tester la précision et l'efficacité de la méthode FETI-DP proposée dans la prédiction de la force optique. La Figure 3a et b présente sa configuration et ses dimensions. Les paramètres constitutifs de l'argent sont toutes des valeurs mesurées prises de [32]. Pour implémenter le schéma FETI-DP, le domaine d'analyse entier est d'abord divisé en 24 sous-domaines. Des maillages plus denses sont requis près de la surface métallique afin de modéliser l'effet d'amélioration du champ local plasmonique. Des éléments tétraédriques sont adoptés pour la discrétisation, ce qui mène à totalement 6,9 × 10⁵ inconnues, incluant 4,1 × 10⁴ inconnues duales et 313 inconnues de coin. La lumière incidente illumine le long de la direction +z, tandis que la direction de polarisation du champ électrique est −x.

Configuration d'une structure nanocapsule d'argent. a Vue 3D. b Vue frontale et dimensions, où R = 30 nm et h = 60 nm.

Premièrement, nous changeons la longueur d'onde de la lumière incidente λ de 200 nm à 400 nm pour simuler les forces optiques exercées sur la nanocapsule. Puisque le FETI-DP fonctionne dans le domaine fréquentiel, les forces optiques sont calculées à 15 points de fréquence d'échantillonnage. La Figure 4 montre la courbe calculée des forces optiques exercées sur la nanocapsule d'argent. Pour indiquer la précision du FETI-DP, les résultats de force optique du FETI-DP sont comparés avec ceux du logiciel commercial Lumerical FDTD Solutions [33], et un bon accord peut être observé.

Résultats des forces optiques exercées sur la nanocapsule d'argent, variant avec la longueur d'onde λ de la lumière incidente, incluant les résultats du FETI-DP et du logiciel commercial FDTD Solutions.

Ensuite, la performance du FETI-DP est testée pour différents nombres de sous-domaines. Nous augmentons le nombre de sous-domaines de 4 à 24 en gardant la densité de discrétisation. Nous assignons chaque processeur pour traiter un sous-domaine. Le Tableau 1 rapporte le temps utilisé pour la construction de l'équation d'interface globale (6) et le temps total de solution. On peut voir que le FETI-DP peut exploiter pleinement les ressources de calcul parallèle et améliorer significativement l'efficacité de la solution. De plus, la précision du FETI-DP avec l'augmentation du nombre de sous-domaines est aussi examinée et rapportée dans le Tableau 1. Ici, la précision est définie par l'erreur relative de norme 2 de la force optique comme δ_OF = ‖OF_i − OF_ref‖/‖OF_ref‖, où OF_i est la force optique utilisant i sous-domaines et OF_ref dénote la force optique de référence utilisant deux sous-domaines. On peut voir que la précision reste presque constante avec l'augmentation du nombre de sous-domaines.

Dimère de nanosphères d'or

Le second exemple analyse le piégeage optique d'une nanosphère d'or en utilisant un dimère de nanosphères d'or. Les effets plasmoniques au gap du dimère peuvent efficacement améliorer la force optique pour le piégeage de nanoparticules. La Figure 5a et b donne la configuration et les dimensions de ce système. Les paramètres constitutifs de l'or sont toutes des valeurs mesurées prises de [32]. Le milieu environnant est l'eau avec un indice de réfraction relatif de n = 1,33. La lumière incidente est une onde plane avec une puissance de 10 mW/µm², la direction de polarisation du champ électrique est +x, et la direction d'incidence est −z. La force optique exercée sur l'objet nanosphère est calculée par la méthode FETI-DP. Pour l'implémentation FETI-DP, le domaine computationnel entier est divisé en 32 sous-domaines et discrétisé par maillages tétraédriques, ce qui résulte en 3,5 × 10⁶ inconnues, incluant 1,6 × 10⁵ inconnues duales et 1738 inconnues de coin.

Configuration d'un système de piégeage optique de dimère de nanosphères d'or dans l'eau. a Vue 3D. b Vue frontale et dimensions, où R = 25 nm, r = 5 nm, et g = 2 nm.

Premièrement, nous testons la performance parallèle du FETI-DP proposé en utilisant divers nombres de processeurs. Le Tableau 2 rapporte le temps de solution pour l'Éq. (6) ainsi que le temps total de solution. De plus, les accélérations pour le calcul parallèle sont aussi fournies dans le Tableau 2. Ici, l'accélération est définie par :

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

où \( {T}_{N_p} \) dénote le temps total d'horloge murale utilisant N_p processeurs. On peut voir que le FETI-DP améliore significativement l'efficacité de la solution et exhibe une bonne accélération parallèle. Pour ce grand nombre d'inconnues, l'utilisation totale de mémoire de tous les processeurs est seulement 57,2 Go.

Ensuite, l'efficacité de l'approche de parcimonie de bas rang est examinée. Avec la parcimonie de bas rang, la matrice de sous-domaine peut être factorisée par algorithme parcimonieux en données et stockée comme matrices parcimonieuses en données. Le temps de construction et l'utilisation de mémoire sont seulement 18 s et 0,5 Go, tandis qu'ils sont 67 s et 1,7 Go par algorithme matriciel conventionnel. On peut voir que nous obtenons 72 % d'économie de temps et 70 % de compression de mémoire. Relié à l'utilisation de mémoire, les PMV subséquents peuvent aussi obtenir 70 % d'économie de temps.

Ensuite, le FETI-DP est testé pour le calcul de la force optique avec la longueur d'onde λ variant de 277 nm à 818 nm. En pratique, les analyses de force optique sous lumière incidente de différentes longueurs d'onde sont souvent nécessaires pour rechercher la longueur d'onde de résonance plasmonique, où l'amélioration drastique du champ se produit et la force optique la plus forte peut être obtenue. Deux cas sont considérés avec la nanosphère située à (0, 0, 20 nm) et (0, 0, −20 nm). La Figure 6a et b trace les forces optiques calculées exercées sur la nanosphère pour différents λ. On peut voir que la force optique maximale se produit à λ = 472 nm, qui est la longueur d'onde de résonance plasmonique. La force optique à cette longueur d'onde de résonance est améliorée de près de 40 fois par rapport à celle à longueur d'onde non-résonance. De plus, la force optique pointe toujours vers le gap du dimère, comme montré en Fig. 6, où l'intensité du champ électrique est la plus forte. C'est aussi la direction de la force de gradient pour piéger l'objet. La Figure 7a et b montre les distributions calculées d'amélioration du champ électrique au plan xoz pour la longueur d'onde non-résonance de λ = 300 nm et la longueur d'onde de résonance de λ = 472 nm, respectivement. On peut voir que l'intensité du champ électrique a été augmentée de presque 250 fois due à l'effet de résonance plasmonique.

Résultats calculés des forces optiques exercées sur la nanosphère dans le système de dimère de nanosphères d'or, variant avec la longueur d'onde λ de la lumière incidente. a La nanosphère objet est située à (0, 0, 20 nm). b La nanosphère objet est située à (0, 0, −20 nm).

Distributions d'amélioration du champ électrique sur le plan xoz pour le système de dimère de nanosphères d'or. a λ = 300 nm (longueur d'onde non-résonance). b λ = 472 nm (longueur d'onde de résonance).

De plus, la force optique et le potentiel optique sont calculés avec la nanosphère se déplaçant de (0, 0, −30 nm) à (0, 0, −17 nm) le long de l'axe z. Puisque le comportement le plus typique et intéressant des forces et potentiels de piégeage sont ceux agissant le long de la direction z, nous considérons ici le potentiel de piégeage axial par intégration le long de l'axe z. Parce que le mouvement de la nanoparticule est restreint à une dimension, la définition d'un potentiel optique est non ambiguë de (16), même si la force optique totale de (15) est non-conservative. Comme montré en Fig. 8a, b, avec la nanosphère se déplaçant vers le gap du dimère, la force optique et la profondeur de potentiel optique augmentent évidemment. À la position de (0, 0, −17 nm), une profondeur de potentiel optique de 4,6 k_BT est produite, ce qui est suffisant pour surmonter le mouvement brownien dans l'eau pour atteindre un piégeage optique stable.

Forces optiques et potentiels optiques exercés sur la nanosphère dans le système de dimère de nanosphères d'or, quand la nanosphère se déplace de (0, 0, −30 nm) à (0, 0, −17 nm). a Les forces optiques. b Les potentiels optiques.

Finalement, nous testons les effets du substrat diélectrique pour cet exemple. Les forces optiques sont calculées avec et sans substrat, respectivement. Pour les deux cas, la nanosphère est située à (0, 0, −20 nm) et la longueur d'onde d'incidence est choisie comme la longueur d'onde de résonance. Pour le cas sans substrat, le résultat calculé de la force optique est |F₀| = 0,769 pN. Pour le cas avec un substrat, le dimère de nanosphères d'or est mis sur un substrat diélectrique avec une épaisseur de 60 nm et une permittivité relative de ε_r = 2,25. Le résultat calculé de la force optique est |F₁| = 0,761 pN. L'erreur relative entre ces deux résultats de forces optiques est d'environ 1,0 × 10⁻², qui est définie comme |F₁ − F₀|/|F₀|.

Dimère de cônes tronqués d'or

Le troisième exemple traite du piégeage optique d'une nanosphère d'or en utilisant un dimère de cônes tronqués d'or. La Figure 9 donne la configuration et les dimensions de ce système. Les paramètres constitutifs de l'or sont pris de [32]. Le substrat diélectrique a une permittivité relative de ε_r = 2,25. Le milieu environnant est l'eau avec un indice de réfraction relatif de n = 1,33. La lumière incidente est une onde plane avec une puissance de 10 mW/µm², la direction de polarisation du champ électrique est +x, et la direction d'incidence est −z. Le domaine computationnel entier est divisé en 32 sous-domaines et discrétisé par maillages tétraédriques, ce qui mène à 3,1 × 10⁶ inconnues, incluant 1,3 × 10⁵ inconnues duales et 1227 inconnues de coin.

Configuration d'un système de piégeage optique de dimère de cônes tronqués d'or basé sur un substrat diélectrique dans l'eau. a Vue 3D. b Vue frontale et dimensions, où UR = 20 nm, LR = 30 nm, h = 35 nm, et g = 2 nm.

Premièrement, nous analysons les forces optiques en changeant λ de 277 nm à 818 nm. La Figure 10 trace les forces optiques calculées exercées sur la nanosphère pour différents λ. La nanosphère est située à (0, 0, 35 nm). On peut voir que la force optique maximale se produit à λ = 464 nm, qui est la longueur d'onde de résonance plasmonique, et la force optique ici est améliorée de près de 30 fois à la longueur d'onde non-résonance. De plus, la force optique totale pointe toujours vers −z, comme montré en Fig. 10, qui est la direction de la force de gradient. Cela confirme que la force de gradient est supérieure à la force de diffusion, qui est une des conditions pour que la nanosphère puisse être piégée stablement. La Figure 11a et b présente les distributions calculées du champ électrique aux longueurs d'onde non-résonance de λ = 300 nm et de résonance de λ = 464 nm, respectivement. On peut voir que l'intensité du champ électrique a été augmentée de presque 500 fois due à la résonance de plasmon de surface localisée.

Résultats calculés des forces optiques exercées sur la nanosphère dans le système de dimère de cônes tronqués d'or, variant avec λ. La nanosphère est située à (0, 0, 35 nm).

Distributions d'amélioration du champ électrique sur le plan xoz pour le système de dimère de cônes tronqués d'or. a λ = 300 nm (longueur d'onde non-résonance). b λ = 464 nm (longueur d'onde de résonance).

Ensuite, la position de la nanosphère est changée à 0, 5 et 35 nm pour observer la force optique. La Figure 12 donne les forces optiques calculées exercées sur la nanosphère, où une composante y évidente de force optique peut être observée, tandis qu'une plus grande composante z de force optique existe. La force optique totale pointe toujours vers la position avec le champ électrique le plus fort pour piéger la nanosphère.

Résultats calculés des forces optiques exercées sur la nanosphère dans le système de dimère de cônes tronqués d'or variant λ. La nanosphère est située à (0, 5 nm, 35 nm).

En outre, nous analysons le potentiel optique avec la nanosphère se déplaçant de (0, 0, 50 nm) à (0, 0, 20 nm) le long de l'axe z. Ici, nous considérons le potentiel de piégeage axial le long de la direction z, qui restreint le mouvement de la nanoparticule à une dimension et mène à une définition non ambiguë du potentiel optique. Les deux forces optiques et potentiels sont calculés. Comme on peut observer de Fig. 13a, b, avec la nanosphère se déplaçant vers le gap du dimère, la force optique et la profondeur de potentiel optique augmentent évidemment. À (0, 0, 20 nm), une profondeur de potentiel optique de 3,8 k_BT est obtenue, ce qui est suffisant pour surmonter le mouvement brownien dans l'eau pour atteindre un piégeage optique stable.

Forces optiques et potentiels optiques exercés sur la nanosphère dans le système de dimère de cônes tronqués d'or, quand la nanosphère se déplace de (0, 0, 50 nm) à (0, 0, 20 nm). a Les forces optiques. b Les potentiels optiques.

Finalement, nous testons les coûts computationnels du FETI-DP en changeant le nombre d'inconnues de 1,0 million à 3,2 millions basé sur différentes tailles de maillage. En pratique, les tests sous différentes densités de maillage sont habituellement nécessaires pour rencontrer différentes exigences de précision. Un tel problème complexe à grande échelle pose de grands défis aux méthodes numériques conventionnelles. Cependant, le FETI-DP peut facilement gérer ce problème. Trente-deux processeurs sont employés pour la simulation FETI-DP, tandis que chaque processeur traite un sous-domaine. Le Tableau 3 rapporte les coûts computationnels du FETI-DP. On peut voir que le FETI-DP exhibe une haute efficacité de simulation et une faible exigence de mémoire.

Conclusion

Une méthode FETI-DP combinée avec parcimonie de bas rang est proposée pour la prédiction et l'analyse du piégeage optique de nanoparticules métalliques. La méthode proposée fournit des problèmes de sous-domaines pleinement découplés, ce qui convertit un problème complexe à grande échelle en une série de problèmes simples à petite échelle. Elle est bien adaptée au calcul parallèle et peut améliorer significativement l'efficacité de la simulation numérique. Les exemples démontrent que la méthode proposée exhibe une excellente performance de calcul à grande échelle et est bien adaptée pour la simulation rapide et précise du piégeage optique à l'échelle nanométrique.

Disponibilité des données et matériaux

Toutes les données générées ou analysées pendant cette étude sont incluses dans cet article.

Abréviations

ABC:

Condition de frontière absorbante

DDL:

Degrés de liberté

FDTD:

Différences finies dans le domaine temporel

MEF:

Méthode des éléments finis

FETI-DP:

Déchirement et interconnexion d'éléments finis dual-primal

TCM:

Tenseur de contrainte de Maxwell

PMV:

Produit matrice-vecteur