Processus d’aimantation dans les nanostructures 53

Champ magnétostatique ou démagnétisant

Ce champ résulte des interactions dipolaires. Contrairement aux précédents, il est à longue portée ; couplant toutes les parties du système, il est sensible à la géométrie globale de l’échantillon. Sa non-localité en fait le champ le plus coûteux à calculer et, de fait, la méthode choisie pour son calcul est déterminante pour l’efficacité du code. Le chapitre 2 (p. 17) est consacré aux méthodes inédites de calcul rapide utilisées dans feeLLGood.
En résumé ici, il s’agit d’un champ irrotationnel [Jac75] dérivant d’un potentiel scalaire ϕ qui vérifie l’équation de Poisson ∆ϕ = ∇ · M . Il présente par ailleurs une discontinuité normale d’amplitude −M · n à la frontière du système. En notant G la fonction de Green pour le laplacien et ∗ la convolution sur Ω, on a alors.

L’ordre deux

L’enjeu est ici de construire un schéma d’ordre deux en temps, permettant une précision équivalente avec des pas de temps beaucoup plus grands. On peut s’inspirer de la méthode de Crank-Nicolson pour l’équation de la chaleur, qui est un θ -schéma avec θ = 1/2 [LP97] ; mais ici, faire θ = 1/2 dans le schéma précédent ne suffit pas. Revenons pour le voir à l’équation (1.30). Pour le moment, seul le laplacien a été décentré ; les autres termes demeurent traités explicitement. En particulier, le dernier terme (|∇m| 2 − H r · m) m, spécifique du micromagnétisme, n’est plus rigoureusement orthogonal au plan tangent K n si l’on y remplace m par m + θ kv. Ainsi deux nouveaux termes d’ordre k apparaissent dans la formulation faible, qui s’écrit : trouver v ∈ K n tel que pour tout w ∈.

Méthodes multipolaires

On va présenter la méthode des multipôles en 2D, ce qui permet d’utiliser la représentation complexe log(z − z ′ ) du potentiel, associant l’affixe z = x + i y au point r ( x , y ). Le développement en série est alors plus simple à obtenir ; en 3D il faut avoir recours auxharmoniques sphériques ce qui n’augmente pas la clarté de l’exposé.

Méthode multipôle hiérarchique

On établit maintenant la structure hiérarchique, une division dyadique du système en cubes : le cube englobant le système est divisé en quatre cubes (en 2D) appelés ses fils dans la terminologie des arbres, et ainsi de suite jusqu’aux cubes terminaux contenant au plus un nombre fixé de sources. On calcule le développement multipolaire de tous les cubes (ie. de la boule circonscrite). Le temps de calcul est M p à chaque niveau, chaque source contribuant aux p coefficients de son cube ; et M p log M au total, en considérant qu’il y a environ log M niveaux.
Le développement d’un cube n’est évidemment pas valable dans les cubes voisins ; il faut aller plus loin. (On entend par « voisins » des cubes de même niveau ayant un point commun.) Deux cubes qui ne sont pas voisins, mais dont les pères le sont sont dits semivoisins. La contribution des semi-voisins est celle calculée par multipôles (voir la Fig. 2.2).
Pour des sources plus proches il faut monter au niveau supérieur (les semi-voisins du fils contenant la cible) ; pour des sources plus éloignées on peut descendre au niveau inférieur.
L’ensemble des semi-voisins de la suite de cubes emboîtés contenant la cible constitue une partition du système privé des voisins du cube terminal. Les sources situées dans cette dernière zone (champ proche) sont calculées à la main, mais elles sont en nombre fixe de sorte que tous les calculs champ proche se font en temps O (N ). En revanche, chaque cube ayant jusqu’à 27 semi-voisins, l’évaluation des développements multipolaires demande 27p opérations par niveau et par cible, soit 27N p log M en tout.
En fin de compte, la méthode multipôle hiérarchique est en O (N log N ), les facteurs limitants étant le calcul et l’évaluation des développements multipolaires, qui nécessitent d’accéder, à chaque niveau, à toutes les sources et toutes les cibles respectivement.

Méthode multipôle rapide (FMM)

La remarque qui va permettre d’accélerer la méthode multipôle hiérarchique en O (N ) est qu’on peut translater un développement du centre d’un cube au centre de son père [GR87, Gre88, CGR88]. Ainsi les développements sont calculés à partir du niveau immédiatement supérieur, en sommant les translatés des quatre fils, sans avoir à remonter jusqu’aux sources. La translation d’un point z 0 à l’origine se fonde sur la relation.

La régularisation intérieure

On se souvient qu’un développement multipolaire n’était valable que suffisamment loin du pôle. Ici, on a un problème analogue : la fonction de Green ne peut être développée telle quelle en série de Fourier à cause de la singularité polaire à l’origine. On va donc effectuer un lissage préliminaire de la fonction de Green sur une boule centrée à l’origine et de rayon ǫ I qu’on appelle la zone de régularisation intérieure. La régularisation est schématisée en 1Dsur la Fig. 2.3. Si l’on note ˜G la régularisée de G, on a donc.

Compromis précision / espace mémoire

On a implémenté dans feeLLGood aussi bien la NFFT que la FMM. Le jeu des préfacteurs fait que la complexité asymptotique, O (N log N ) ou O (N ), ne préjuge pas de la méthode qui sera plus efficace sur les maillages typiques de 10 4 à 10 6 nœuds. Le paramètre qui commande la précision en NFFT est le nombre de modes (en FMM, c’est le nombre de termes des développements).
Il faut remarquer qu’une augmententation du nombre de modes n’accroît pas seulement le temps calcul, mais également l’espace mémoire utilisé. Chacune des deux NFFT correspondant à la double somme (2.41) fait usage, pour l’étape FFT, de deux plans FFTW (un en entrée, un en sortie). Pour éviter tout problème de repliement du spectre, les FFT sont effectuées avec un facteur 2 de sur-échantillonnage, dans chaque dimension, par rapport aux modes demandés par l’utilisateur.
On doit donc réserver assez de mémoire pour 4 × 8n3 nombres complexes, auxquels s’ajoutent ν N réels pour la correction de champ proche selon (2.35), et (2m + 2)3 M réels pour la diffusion sur le réseau. En règle générale, on a observé que la méthode NFFT demande 2 à 4 fois plus d’espace mémoire que la FMM pour un calcul équivalent.

Optimisation pour les systèmes anisotropes

On se souvient que dans la méthode NFFT, la fonction de Green régularisée est échantillonnée sur un cube contenant le système. Pour un système plat, ie. de rapport d’aspect élevé, cela signifie que la grande majorité des points d’échantillonnage se trouve en dehors du système. Or le comportement de la fonction de Green dans le vide n’est pas d’un grand intérêt. Pour optimiser la méthode précédente, l’idée est d’échantillonner G sur un rectangle adapté aux dimensions du système, de rapport d’aspect voisin (voir la Fig. 2.5).
Pour faire cette optimisation, il faut néanmoins renoncer à l’une des deux régularisations vues précédemment : la régularisation de bord. En effet, celle-ci s’effectue à l’extérieur d’une boule contenant le système et contenue dans la cellule, cubique en NFFT, de périodisation. L’équivalent anisotrope de cette boule est un ellipsoïde ; cependant on a vu que les problèmes de raccordement des fonctions radiales ne peuvent être résolus que sur une sphère.
L’absence de régularisation de bord entraîne une discontinuité de la dérivée première de la fonction de Green périodisée, aux points de recollement des périodes. Il est toujours possible dans la version anisotrope de choisir une continuité C p pour le raccordement intérieur ; en revanche le raccordement au bord du rectangle est forcément C 0 . Comment la classe du raccordement influe-t-elle sur la précision du calcul NFFT ?
On a fait des expériences numériques comparant le calcul NFFT avec différentes classes de raccordement au calcul direct (quadratique). La notation C p − C q signifie un raccordement intérieur C p et un raccordement de bord C q . On a choisi un système à aimantation non uniforme, qui présente aussi bien des charges volumiques que surfaciques : la sphère de Halbach. Il s’agit en réalité d’une couronne sphérique dont l’aimantation, planaire, est donnée par m = cos θ u r +sin θ u θ en coordonnées polaires. On peut se représenter l’aimantation dans une section suivant son plan, comme tournant deux fois plus vite que l’angle polaire (voir la Fig. 2.6).
Notons qu’on dispose pour ce système d’une expression analytique du potentiel (le champ est même uniforme dans la cavité), mais on n’en fait pas usage ici de façon à ne pas confondre l’erreur due au maillage et l’erreur NFFT. On prend donc pour référence, sur un maillage donné, le potentiel obtenu par la sommation directe de (2.20).

Problème standard no4 du NIST

Présentation du problème

On s’intéresse maintenant à un système à distribution d’aimantation non uniforme. Il s’agit d’une mince plaquette de Permalloy sous champ planaire, de dimensions L = 500 nm, l = 125 nm et t = 3 nm (voir la Fig. 3.9). L’état initial de la simulation est l’état S (voir la Fig. 3.13) obtenu en diminuant progressivement un champ de direction (111) depuis la saturation jusqu’à la rémanence. Les paramètres matériau sont A = 1.3 × 10−11 J/m et M s = 8.0 × 10 5 A/m. On prend α = 0.02 ; on s’attend à une dynamique fortement non-linéaire.
Deux champs extérieurs sont envisagés dans le plan ( x , y ) de la plaquette, l’axe des x étant celui de la longueur :
– le champ 1 d’angle azimutal φ = 170 ◦ et d’amplitude 25 mT ; – le champ 2 d’angle azimutal φ = 190 ◦et d’amplitude 36 mT.
Le champ extérieur va sortir l’aimantation du plan et créer des charges de surface, le champ démagnétisant induit entraînant la rotation du système. A partir de l’état S, la première composante du champ extérieur fait tourner les bords gauche et droit du système et la deuxième composante en fait tourner le cœur. Le cœur tourne donc dans des sens opposés suivant qu’on applique le champ 1 ou le champ 2. Avec le champ 1, il tourne dans le même sens que les bords ; et dans l’autre sens avec le champ 2. C’est pourquoi leproblème 2 est sensiblement plus compliqué.
Dans la suite on compare sur différents maillages, les résultats obtenus par feeLLGood et le code différences finies GL-FFT [Kev98]. D’autres codes différences finies sont compa-44 FIG. 3.10: Comaparaison de codes différences finies sur le site du NIST : évolution de 〈m y 〉 sous champ 1 (g.) et 2 (d.). rés sur le site du NIST (voir la Fig. 3.10).
En revanche, rares sont les résultats éléments finis publiés sur la question. Dans [TSSF01] on pourra voir l’évolution en énergie pour le champ 1 mais on la cherchera en vain pour le champ 2. Ce problème plus difficile est ici analysé en détail.
La référence utilisée est le code GL-FFT sur un maillage en cellules de 1 × 1 × 3 nm.
Le système ayant un rapport d’aspect de l’ordre de 100, on utilise la méthode des NFFT anisotropes pour les calculs de champ démagnétisant. C’est d’ailleurs un cas typique où les codes utilisant FEM-BEM sont peu efficaces.

Comparaison des champs 1 et 2 et maillages tétraédriques

On commence par tester le schéma d’ordre un. On sait, grâce à l’expérience du sphéroïde macrospin, que compte tenu de la valeur de α le critère géométrique avec d = 0.01 est insuffisant. Au lieu de diminuer d , on réalise ici une autre possibilité : prendre un critère énergétique plus sévère. En l’espèce la dissipation numérique est bornée au dixième de la dissipation physique.
Le système est beaucoup plus gros que le sphéroïde : les dimensions L et l sont de l’ordre de la centaine de longueurs d’échange. On s’attend donc à devoir utiliser des maillages plus lourds. Pour déterminer le maillage adapté, on a effectué pour chaque champ appliqué la simulation sur deux maillages de respectivement 86 000 et 168 000 tétraèdres.
On constate que la dynamique sous champ 1 est déjà bien décrite avec le maillage à 86 000 tétraèdres (Fig. 3.11) ; ce n’est pas le cas en revanche de la dynamique sous champ 2, qui s’écarte sensiblement de la référence vers t = 0.5 ns (Fig. 3.12). Cet instant correspond à l’apparition de zones d’aimantation perpendiculaire (configuration D de la figure 3.14).

Optimisations : maillages extrudés et schéma d’ordre deux

Le maillage à 168 000 tétraèdres donne de bons résultats, y compris sous champ 2 ; les configurations conservent une symétrie remarquable par rapport à l’origine, jusqu’après la nucléation (Fig.3.14). Cependant, un maillage de cette taille demande pas moins de 16 s pour l’exécution d’une itération. On se demande s’il est possible de décrire le système avec moins d’éléments, en choisissant par exemple une forme d’éléments plus adaptée au problème.
On a donc généralisé feeLLGood de façon à exploiter plusieurs types de maillages : non seulement tetraédriques, mais aussi pentaédriques et hexaédriques (voir Fig. 3.15).
L’avantage de ces maillages est qu’ils peuvent s’obtenir par extrusion, et exhiber ainsi une invariance par translation suivant une coordonnée. Pour un système plat, il devient possible d’introduire une ou plusieurs couches en épaisseur du maillage, ce qui peut être décisif.
En effet, l’erreur sur le champ démagnétisant due à l’intégration numérique, se concentre typiquement aux nœuds de surface.

Sélection de chiralité dans un plot monovortex

Les études de Raabe et al. sur des plots de cobalt ainsi que Shinjo et al. [SOH + 00, RPS + 00] sur des plots de Permalloy ont mis en évidence, dans la gamme de diamètres 0.1 à 1 µm et à la rémanence, l’état dit vortex. Il s’agit d’un état à fermeture de flux, favorable du point de vue démagnétisant : l’aimantation circule le long des bords du système, et est radialement invariante excepté dans une zone centrale de diamètre 10 à 20 nm appelée le cœur de vortex (voir la Fig. 4.1 à g.). Si la circulation se fait dans le sens horaire (resp. anti-horaire), on parle de chiralité droite (resp. gauche). Dans le cœur de vortex, l’aimantation sort du plan pour éviter l’occurrence d’un point singulier et respecter ainsi la contrainte de norme uniforme. Un axe nord-sud normal au plan étant défini, on parle de polarité nord ou sud suivant l’orientation de l’aimantation dans cette zone.
Polarité et chiralité sont deux degrés de liberté qu’il est possible de manipuler indépendamment [ZZP00, KOKS01], permettant en théorie de stocker deux bits d’informations dans un plot monovortex. Le retournement de polarité, par l’application d’un champ hors du plan, est relativement bien compris [OSO + 02] ; en revanche les processus qui commandent la chiralité dépendent de la forme du plot et sont expliqués jusqu’à présent defaçon heuristique [SHZ01].

Cas du plot circulaire tronqué

Les travaux de Schneider et al. [SHZ00] ont montré qu’il n’est pas possible de commander la chiralité dans un plot parfaitement circulaire ; il est nécessaire pour cela d’introduire une anisotropie géométrique. Schneider a ainsi étudié des plots de Permalloy à section tronquée suivant une corde, exhibant donc une surface courbe et une facette droite [SHZ01].
En saturant le système par un champ extérieur parallèle à la facette et en laissant décroître le champ, on arrive à rémanence à l’état vortex. La microscopie de Lorentz montre alors que l’aimantation au niveau de la facette est toujours antiparallèle à l’état saturé initial. Autrement dit, la chiralité est sélectionnée par le signe du champ appliqué.
L’explication heuristique de Schneider est la suivante. À mesure que le champ diminue, les lignes d’aimantation s’incurvent pour suivre la surface du plot. L’aimantation est piégée en premier lieu par la surface courbe, la plus longue ; le retournement s’opère au voisinage de la facette. On peut donc supposer qu’un vortex entre par la facette et finit par se stabiliser à l’intérieur du système.

Dépendance angulaire dans un plot facetté de symétrie C3

On veut à présent traiter le cas de champs appliqués non parallèles à la facette, et étudier la dépendance angulaire de la chiralité. On s’attend à ce que la position du champ relativement aux axes de symétrie du système joue un rôle important. Le système choisi est un plot hexagonal de cobalt à trois sur-facettes conférant une symétrie C 3 (trigonale).
Le diamètre du plot est 250 nm et la hauteur 50 nm. Paramètres matériau : A = 3 × 10 −11 J/m et M = 1.42 × 10 6A/m.
La Fig. 4.3 montre le processus de formation du vortex pour un champ appliqué d’angle polaire θ = 25
◦ avec la normale au plan du plot et angle azimutal φ = 10
◦ avec la facette supérieure.
C’est un vortex de polarité nord et chiralité gauche qui entre par la facette. La polarité sud est défavorisée énergétiquement par la composante hors du plan du champ appliqué ; on arrive ainsi, à la rémanence, à une chiralité gauche (Fig. 4.3).
La Fig. 4.4 montre le processus symétrique, avec un champ appliqué suivant φ = 170◦.
La nucléation se développe toujours à partir de la facette supérieure, mais le vortex favorisé, de polarité nord, a cette fois une chiralité droite.
On a effectué des simulations pour une variété d’angles φ avec un pas de 10 ◦ , avec toujours α = 1. Sur la Fig. 4.5, les chiralités obtenues sont comparées aux chiralités mesurées par Massebœuf et al. [MFT+ 10] en microscopie de Lorentz.
On observe un très bon accord des chiralités prévues par la simulation avec celles observées expérimentalement. Les simulations permettent de distinguer clairement six secteurs de chiralités alternativement gauche et droite. Ce résultat reflète la division du système en six régions par la symétrie C 3 (trois axes). Il est également visible par microscopie, si l’on excepte quelques points de désaccord facilement imputables aux imperfections des plots expérimentaux.
Les mêmes simulations ont été effectuées en variant cette fois l’angle polaire. La Fig. 4.6compare les chiralités obtenues respectivement pour θ = 25
◦ et θ = 90
◦ c’est-à-dire un champ planaire.
On observe un léger décalage des chiralités dans le sens horaire pour θ = 90 ◦ . Néanmoins, on voit en comparant les Fig. 4.7 et 4.3 que le processus de nucléation, avec compétition de deux vortex, est fondamentalement inchangé.

Table des matières

1 Micromagnétisme et méthodes numériques 
1.1 Genèse du micromagnétisme
1.2 L’équation de Landau-Lifchitz-Gilbert
1.3 Champs et énergies en jeu
1.3.1 Champ d’échange
1.3.2 Champ d’anisotropie magnétocristalline
1.3.3 Champ magnétostatique ou démagnétisant
1.3.4 Champ extérieur ou Zeeman
1.4 Méthodes numériques
1.4.1 Méthode des différences finies
1.4.2 Méthode des éléments finis
1.5 Formulations faibles et schémas en temps
1.5.1 Schéma explicite et renormalisation
1.5.2 Le θ -schéma
1.5.3 L’ordre deux
2 Calcul du champ magnétostatique 
2.1 Les équations de la magnétostatique
2.2 Méthode des éléments finis (FEM)
2.3 Méthode mixte élements finis – intégrale de frontière (FEM-BEM)
2.4 Généralités sur les méthodes de Green
2.5 Méthodes multipolaires
2.5.1 Méthode multipôle hiérarchique
2.5.2 Méthode multipôle rapide (FMM)
2.6 Méthode de Fourier non-uniforme (NFFT)
2.6.1 La régularisation intérieure
2.6.2 La régularisation de bord
2.6.3 La transformée de Fourier non uniforme
2.6.4 Compromis précision / espace mémoire
2.6.5 Optimisation pour les systèmes anisotropes
3 Validation du code – Application à des problèmes standard 
3.1 Relaxation d’un sphéroïde macrospin
3.2 Problème standard no4 du NIST
3.2.1 Présentation du problème
3.2.2 Comparaison des champs 1 et 2 et maillages tétraédriques
3.2.3 Optimisations : maillages extrudés et schéma d’ordre deux
4 Processus d’aimantation dans les nanostructures 53
4.1 Sélection de chiralité dans un plot monovortex
4.1.1 Cas du plot circulaire tronqué
4.1.2 Dépendance angulaire dans un plot facetté de symétrie C3
4.1.3 Nucléation en faible amortissement
4.2 Résonance ferromagnétique du plot monovortex
4.3 Hystérésis de Néel caps dans un plot de fer
5 Perspectives : étude d’un oscillateur hyperfréquence 
5.1 Contexte
5.2 Résultats
A Le théorème de Bartels
Bibliographie

Lire le rapport complet