# Calcul par éléments finis de diffusion élastique d'ondes électroniques : problèmes direct et inverse

Denis Aubry, Ann-Lenaïg Hamon, Guillaume Puel

► To cite this version:

Denis Aubry, Ann-Lenaïg Hamon, Guillaume Puel. Calcul par éléments finis de diffusion élastique d'ondes électroniques: problèmes direct et inverse. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01418157

HAL Id: hal-01418157

https://hal.science/hal-01418157v1

Submitted on 16 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Calcul par éléments finis de diffusion élastique d'ondes
électroniques : problèmes direct et inverse

D. Aubry, A.-L. Hamon, G. Puel

Laboratoire MSSMat (Ecole Centrale Paris/CNRS UMR8579)
Grande Voie des Vignes, 92290 Châtenay-Malabry
{denis.aubry,ann-lenaig.hamon,guillaume.puel}@ecp.fr

**Résumé** – Nous nous intéressons ici au cas de l'interaction élastique électron-matière survenant dans un microscope électronique en transmission. À partir du potentiel caractérisant l'échantillon observé, on montre que le problème direct permettant d'obtenir la fonction d'onde électronique à la sortie de l'échantillon peut s'écrire comme une équation de Helmholtz résolue par la Méthode des Eléments Finis. Une formulation du problème inverse qui a pour but de retrouver le potentiel de l'échantillon à partir de la fonction d'onde mesurée est également présentée.

**Mots clés** – diffusion élastique, MET, problème inverse, état adjoint, régularisation itérative

# 1 Introduction

Les problèmes de diffusion d'ondes dans un milieu hétérogène interviennent dans de très nombreuses applications et à différentes échelles (sismique, ultrasons, ...) et sont de plus en plus souvent résolus numériquement par la Méthode des Eléments Finis (MEF).

Nous nous intéressons ici au cas de l'interaction élastique électron-matière telle qu'on peut l'observer dans un microscope électronique en transmission (MET)[9]. L'objectif associé est de déterminer la structure d'un échantillon de matière en observant le résultat de la diffusion d'une onde électronique par l'échantillon.

Dans le problème direct, on suppose que le potentiel créé par l'échantillon est connu à l'échelle électronique, et on calcule à l'échelle atomique la fonction d'onde électronique obtenue par diffusion élastique à la sortie de l'échantillon. En introduisant différentes hypothèses simplificatrices mais non restrictives, ce calcul peut être accompli avec un coût raisonnable dans un cas bidimensionnel à l'aide de la MEF. Il est ensuite possible de confronter le résultat obtenu numériquement avec les informations de diffusion obtenues expérimentalement dans le MET.

Cette confrontation peut permettre, par problème inverse, de déterminer le potentiel de l'échantillon supposé inconnu à partir de la fonction d'onde mesurée expérimentalement. Or les solutions d'un problème inverse sont généralement instables et ne sont pas uniques ; nous avons montré précédemment dans [8] que la recherche d'un champ de propriétés variable spatialement était difficile à régulariser, et en particulier que des régularisations de type Tikhonov s'avéraient insuffisantes. C'est pourquoi nous proposons ici une stratégie analogue à celle de [2] et qui repose sur une discrétisation du champ recherché indépendante du maillage utilisé pour le calcul de diffusion électronique. Cette discrétisation, initialement grossière pour assurer une bonne régularisation du problème, est ensuite affinée de manière itérative par l'utilisation d'estimateurs d'erreur classiquement utilisés en adaptation de maillage.
## 2 Cadre théorique du problème direct de diffusion électronique

Dans notre cas, qui correspond à une condition expérimentale classique, l'électron incident est représenté par une fonction d'onde $\psi_i(\mathbf{x}) = \tilde{\psi}_i \exp(i\mathbf{k}_i \cdot \mathbf{x})$ représentée comme une onde plane associée à un vecteur d'onde $\mathbf{k}_i$ et vérifiant dans un domaine vide $\Omega_e$ l'équation de Schrödinger suivante exprimée en unités atomiques (distances en Bohr et énergies en Hartree) :

$$-\frac{1}{2}\Delta\psi_i = E_i\psi_i$$

où $E_i = ||\mathbf{k}_i||^2/2$ est l'énergie de l'électron incident et $\Delta$ l'opérateur laplacien.

Le domaine de l'échantillon, noté $\Omega_s$, est caractérisé par un potentiel $V_s$ associé aux électrons et noyaux de l'échantillon, qui décroît rapidement en dehors de $\Omega_s$.

<figure><img src="image_2.png" /><figcaption>Figure 1 - Domaine considéré</figcaption></figure>

Le système complet à considérer doit faire intervenir à la fois l'électron incident et les particules de l'échantillon au sein de l'équation de Schrödinger générale suivante :

$$-\frac{1}{2}(\Delta_{\mathbf{x}} + \Delta_{\mathbf{x}_s})\psi_{es}(\mathbf{x}, \mathbf{x}_s) + (V_{es}(\mathbf{x}, \mathbf{x}_s) + V_s(\mathbf{x}_s))\psi_{es}(\mathbf{x}, \mathbf{x}_s) = E\psi_{es}(\mathbf{x}, \mathbf{x}_s)$$

où $\psi_{es}$ est la fonction d'onde totale du système d'énergie totale $E$, et le potentiel est décomposé comme la somme de potentiels $V_s$ et $V_{es}$ exprimant l'interaction de l'échantillon avec lui-même et avec l'électron incident respectivement ; $\mathbf{x}$ et $\mathbf{x}_s$ correspondent à la description de la position de l'électron incident et de toutes les autres particules respectivement.

Pour découpler les termes liés à l'électron incident des termes liés à l'échantillon, on peut montrer [11] que l'on peut approximer la fonction d'onde totale $\psi_{es}$ comme le produit de deux fonctions d'ondes $\psi_s$ et $\psi_e$ associées à l'échantillon et l'électron incident respectivement, de sorte que $\psi_e$ doit vérifier l'équation suivante :

$$-\frac{1}{2}\Delta_{\mathbf{x}}\psi_e(\mathbf{x}) + V(\mathbf{x})\psi_e(\mathbf{x}) = (E - E_s)\psi_e(\mathbf{x})$$

où $E_s$ est l'énergie de l'échantillon et $V(\mathbf{x}) = \int_{\Omega_s} V_{es}(\mathbf{x}, \mathbf{x}_s) ||\psi_s(\mathbf{x}_s)||^2 d\mathbf{x}_s$.

Si on décompose la fonction d'onde électronique $\psi_e$ comme la somme de l'onde incidente $\psi_i$ connue et de l'onde transmise $\psi_d$ à déterminer, et si on prend en compte les ordres de grandeur des énergies en présence ($||V|| \ll E_i \approx E - E_s$), on aboutit à l'équation de Helmholtz suivante :

$$\frac{1}{2}\Delta\psi_d + E_i\psi_d = V\psi_i \quad (1)$$

Dans le cas où on tient pas compte des défauts des lentilles du microscope, c'est (à un grandissement et une rotation près) le carré du module de la fonction d'onde électronique, $||\psi_e||^2 = ||\psi_i + \psi_d||^2$, que l'on observe sur un plan virtuel $\Sigma_m$ situé à la sortie de l'échantillon $\Omega_s$ comme représenté sur la Figure 1.
# 3 Résolution numérique du problème direct de diffusion électronique

## 3.1 Approche classique vs. calcul Eléments Finis

L'approche classique consiste à considérer la fonction de Green $g$ associé au problème d'Helmholtz (1) à résoudre :

$$g(\mathbf{r}) = \frac{\exp(i\mathbf{k} \cdot \mathbf{r})}{r}$$

Après multiplication de l'équation (1) par $g$ et intégration par parties, on aboutit à l'expression suivante pour l'onde transmise :

$$\psi_d(\mathbf{x}) = \int_{\mathbb{R}^3} V(\mathbf{y})\psi_i(\mathbf{y})g(\mathbf{x} - \mathbf{y})d\mathbf{y} = ((V\psi_i) * g)(\mathbf{x})$$

Dans l'expression précédente, les réflexions qui vont et viennent entre les atomes ne sont alors pas prises en compte.

La résolution se fait alors de façon approchée par l'application d'une méthode « multi-tranches » [5][6][10][12] qui subdivise l'échantillon en plusieurs tranches fonctions de l'épaisseur, la taille de ces tranches étant choisie de façon plus ou moins empirique.

Plutôt que d'utiliser l'approche précédente, on désire résoudre le plus fidèlement possible l'équation (1) en considérant l'ensemble des ondes interagissant, sans faire les hypothèses de l'approche classique. On souhaite aussi pouvoir étudier le cas d'échantillons non périodiques, car présentant un défaut localisé par exemple.

Pour cela, on souhaite utiliser la MEF. Une première difficulté est liée au fait que le problème initial (1) est posé sur un domaine $\Omega_e$ de dimensions infinies. Pour contrer ce problème, on peut utiliser le fait que l'onde recherchée est, loin de l'échantillon, de type évanescent, ce qui nous permet de borner le domaine $\Omega_e$ par une frontière, notée $\Sigma_\infty$ sur la Figure 1, sur laquelle est imposée la condition évanescente suivante [7] :

$$\frac{\partial \psi_d}{\partial \mathbf{n}} = i \|\mathbf{k}_i\| \psi_d \qquad (2)$$

où $\mathbf{n}$ désigne le vecteur normal unitaire sortant le long de la frontière $\Sigma_\infty$.

Malgré ceci, la principale difficulté restante reste le caractère haute fréquence du problème à résoudre. Ainsi, l'ordre de grandeur de la longueur d'onde incidente est $\lambda_i = 0,05$ u.a. tandis que l'épaisseur totale de l'échantillon est de l'ordre de $5000$ u.a. (pour des distances interatomiques d'environ $5$ u.a.) ; en faisant l'hypothèse de $10$ degrés de liberté EF par longueur d'onde, on arrive à un total d'environ $10^9$ degrés de liberté pour ne mailler qu'une seule cellule cristalline, ce qui serait rédhibitoire sans traitement préalable tel que celui proposé ci-dessous.

## 3.2 Approximation paraxiale

Au lieu de chercher la fonction d'onde $\psi_d$, l'approximation paraxiale consiste à chercher l'inconnue $\tilde{\psi}_d$ telle que :

$$\psi_d(\mathbf{x}) = \tilde{\psi}_d(\mathbf{x}) \exp(i\mathbf{k}_i \cdot \mathbf{x}) \qquad (3)$$

où $\mathbf{k}_i$ est le vecteur d'onde incident, dirigé selon l'axe du microscope dans le cas d'une illumination parallèle. Même si, formellement, cette approximation n'implique pas que le champ diffracté est de direction proche de celle du champ incident, cette condition est compatible avec les conditions expérimentales associées au MET, dans lequel tous les rayons diffractés avec un angle supérieur à $30$ mrad environ sont tronqués par la fonction de transfert du microscope.
En introduisant l'approximation (3) dans l'équation (1), on obtient :

$$ \frac{1}{2} \Delta \tilde{\psi}_d + i \mathbf{k}_i \cdot \nabla \tilde{\psi}_d = V \tilde{\psi}_i \text{ dans } \Omega_e \quad (4) $$

car $E_i = ||\mathbf{k}_i||^2/2$. De même, la condition d'évanescence (2) devient, avec (3) :

$$ \frac{\partial \tilde{\psi}_d}{\partial \mathbf{n}} = i(||\mathbf{k}_i|| - \mathbf{k}_i \cdot \mathbf{n}) \tilde{\psi}_d \text{ sur } \Sigma_{\infty} \quad (5) $$

# 4 Formulation du problème inverse de diffusion électronique

Le problème inverse consiste à déterminer dans $\Omega_e$ le champ spatial de potentiel $V$ tel que l'intensité du champ électronique $||\tilde{\psi}_i + \tilde{\psi}_d||^2$ calculée numériquement à partir de $V$ sur le plan virtuel $\Sigma_m$ corresponde au mieux à l'intensité $||\psi_m||^2$ mesurée à la sortie de l'échantillon. On espère ainsi pouvoir localiser un défaut dans l'échantillon étudié.

Pour cela, on introduit la fonction coût suivante, proposée dans [3][4] :

$$ J(V) = \frac{1}{4} \int_{\Sigma_m} \left( ||\tilde{\psi}_i + \tilde{\psi}_d||^2 - ||\psi_m||^2 \right)^2 + \frac{\alpha}{2} \int_{\Omega_v} (V - V_0)^2 \quad (6) $$

où $\alpha$ est un paramètre de régularisation à fixer, et $V_0$ un potentiel choisi a priori, supposé proche du champ recherché : par exemple, le potentiel associé au cristal parfait dans le cas de la localisation d'un défaut. Dans ce cas, le problème inverse consiste à trouver le défaut grâce à la modification associée du potentiel.

## 4.1 Minimisation de la fonction coût par la méthode de l'état adjoint

La résolution du problème inverse se fait par l'introduction d'un problème d'état adjoint, proche du problème direct, qui permettra d'exprimer directement le gradient de la fonction coût (6). L'état adjoint $z$ vérifie les équations suivantes :

$$ \frac{1}{2} \Delta z + i \mathbf{k}_i^* \cdot \nabla z = 0 \text{ dans } \Omega_e \quad (7) $$

$$ \frac{\partial z}{\partial \mathbf{n}} = -i(||\mathbf{k}_i|| + \mathbf{k}_i^* \cdot \mathbf{n})z \text{ sur } \Sigma_{\infty} \quad (8) $$

$$ \frac{1}{2} \left[ \left[ \frac{\partial z}{\partial \mathbf{n}} \right] \right] = \left( ||\tilde{\psi}_i + \tilde{\psi}_d||^2 - ||\psi_m||^2 \right) (\tilde{\psi}_i + \tilde{\psi}_d) \text{ sur } \Sigma_m \quad (9) $$

où $*$ et [[·]] désignent le conjugué et le saut respectivement. L'état adjoint s'interprète comme un problème de rétropropagation d'onde.

Pour minimiser la fonction coût, la variation de (6) est calculée :

$$ \delta J = \int_{\Sigma_m} \left( ||\tilde{\psi}_i + \tilde{\psi}_d||^2 - ||\psi_m||^2 \right) \operatorname{Re} \left( (\tilde{\psi}_i + \tilde{\psi}_d) \delta \tilde{\psi}_d^* \right) + \alpha \int_{\Omega_v} (V - V_0) \delta V $$

et on peut montrer qu'en prenant la variation de (4) et après plusieurs intégrations par parties faisant intervenir (5)-(7)-(8)-(9), la variation de la fonction coût s'exprime comme :

$$ \delta J = \int_{\Omega_v} - \operatorname{Re}(\tilde{\psi}_i z^*) \delta V + \alpha \int_{\Omega_v} (V - V_0) \delta V $$

et comme la détermination du minimum de la fonction coût est liée à $\delta J = 0 \forall \delta V$, on a :

$$ \operatorname{Re}(\tilde{\psi}_i z^*) = \alpha(V - V_0) \text{ dans } \Omega_e \quad (10) $$
Ainsi, la résolution du problème inverse consiste à résoudre le système couplé d'équations (4)-(7)-(10) d'inconnues {$\tilde{\psi}_d, z, V$} associé aux conditions limites (5)-(8)-(9).

## 4.2 Stratégie itérative à maillages adaptatifs

De même que nous avons choisi une approche par Eléments Finis pour résoudre le problème direct, le système (4)-(7)-(10) est résolu par discrétisation des équations à l'aide d'une base Eléments Finis. Toutefois, le potentiel dont l'identification est associée à la résolution du système décrit précédemment est un champ scalaire variant spatialement : le déterminer revient à déterminer les valeurs de ce potentiel à chaque noeud du maillage utilisé, ce qui rend le problème particulièrement mal posé si le maillage présente un nombre conséquent de degrés de liberté.

C'est pourquoi nous proposons de suivre la stratégie utilisée dans [8] où la discrétisation spatiale du champ à identifier est réalisée sur un maillage différent du maillage permettant la résolution des problèmes direct et adjoint (4)-(7). En effet, de même que Bangerth et al dans [2], nous avons pu constater dans [8] que l'utilisation d'un maillage grossier pour discrétiser le potentiel recherché et ainsi résoudre l'équation de compatibilité (10) agit comme une régularisation supplémentaire du problème inverse, car le potentiel à identifier est recherché dans un espace plus petit qui tend à restreindre l'étendue du potentiel.

Cette régularisation peut alors s'effectuer de manière itérative, dès lors qu'on décide de raffiner progressivement le maillage utilisé pour la discrétisation du potentiel recherché, de façon à améliorer la qualité du champ identifié ; ceci est possible grâce à l'utilisation de méthodes classiques d'adaptation de maillage, où les éléments qui contribuent le plus à l'erreur caractérisée par la fonction coût (6) sont raffinés de manière automatique.

# 5 Exemple numérique

L'exemple choisi à titre d'illustration est un échantillon bidimensionnel de fer-□ pur. La Figure 2 montre la simulation directe obtenue avec 600 000 ddl dans un MET de 200 keV avec $||k_i|| = 130$ rad/(u.a.), où l'échantillon d'épaisseur 2 nm est orienté selon la direction [001], et où le potentiel $V$ associé au cristal est simulé à l'aide du potentiel de Yukawa [1] (superposé sur la Figure 2 de façon à visualiser la position des noyaux de l'échantillon). Comme l'équation (4) est linéaire en $\tilde{\psi}_i$, amplitude de l'onde plane incidente, et en $V_a$, amplitude du potentiel choisi, ces deux valeurs ont été fixées à 1. De plus, un coefficient d'absorption de valeur $5.10^{-3}$ a été introduit pour améliorer la convergence du calcul.

<figure><img src="image_3.png" /><figcaption>Figure 2- Calcul direct sur un échantillon de fer-□</figcaption></figure>
# 6 Conclusion

Nous avons montré la faisabilité de la résolution du problème direct de diffusion électronique à l'aide de la MEF. Ceci nous permet aussi de simuler des données expérimentales avec lesquelles nous avons mené plusieurs premiers essais de résolution du problème inverse.

<figure><img src="image_4.png" /><figcaption>Figure 3- (gauche) Donnée expérimentale simulée (cristal avec lacune) (droite) Identification associée de l'écart au potentiel du cristal parfait</figcaption></figure>

# Références

[1] N. W. Ashcroft, N. D. Mermin. *Physique des Solides*, EDP Sciences (Les Ulis), 2002.

[2] W. Bangerth, A. Joshi, E. MSevick-Muraca. Adaptive FE methods for increased resolution in fluorescence optical tomography. *Progr. Biomed. Optics Imag.* **6**, pp. 318-329, 2007.

[3] L. Beilina, C. Clason. An adaptive hybrid FEM/FDM method for an inverse scattering problem in scanning acoustic microscopy. *J. SIAM Sci. Comp.* **28**:1, pp. 382-402, 2006.

[4] L. Beilina, C. Johnson. A posteriori error estimation in computational inverse scattering. *J. Mathematical Models and Methods in Applied Sciences* **15**:1, pp. 23-37, 2005.

[5] K. Ishizuka, N. Uyeda. A new theoretical and practical approach to the multislice method. *Acta Cryst.* **A 33**, pp. 740-749, 1977.

[6] E. J. Kirkland. *Advanced Computing in Electron Microscopy*, Plenum (New York), 1998.

[7] A. V. Popov. Computation of paraxial wave fields using transparent boundary conditions. *Computational Mathematics and Mathematical Physics* **46**:9, pp.1595-1600, 2006.

[8] G. Puel, D. Aubry. *Identification of a spatial field of material properties with adaptive regularization and meshes*. Proceedings of IDETC/CIE 2008, ASME, CD-ROM, 2008.

[9] A. R. Smith, L. Eyring. Calculation, display and comparison of elecron microscope images modeled and observed. *Ultramicroscopy* **8**, pp. 65-78, 1982.

[10] P. A. Stadelmann. EMS – a software package for elecron diffraction analysis and HREM image simulation in materials science. *Ultramicroscopy* **21**, pp. 131-146, 1987.

[11] Z. L. Wang. *Elastic and inelastic scattering in electron diffraction and imaging*, Plenum (New York), 1995.

[12] D. B. Williams, C. B. Carter. *Transmission electron microscopy*, Plenum (New York), 1996.