Structure électronique des interfaces Co(OOOl)/MoS2 et Ni(lll)/WSe2 pour l’injection de spin dans un semi-conducteur bidimensionnel

Les différentes fonctionnelles 

La fonctionnelle de la densité locale (LDA) L’étape cruciale d’un calcul DFT est de trouver la fonctionnelle d’échange et de corrélation, pour laquelle on doit utiliser des approximations qui donnent des résultats acceptables. Si celle-ci était connue, alors la seule difficulté serait de résoudre numériquement l’équation de Kohn et Sham. Or, la fonctionnelle d’échange et de corrélation universelle n’est pas connue de façon exacte, il nous faut donc l’approximer. La première des fonctionnelles utilisées fut l’approximation de la densité locale LDA (Local Density Approximation), proposée par Kohn et Sham en 1965 [16]: E LD xc = Z ρ(r)xc(ρ(r))dr (1.14) où xc(ρ) est l’énergie d’échange et de corrélation par électron pour un gaz d’électrons uniforme de densité ρ. La fonctionnelle d’échange et de corrélation peut être partitionnée en deux parties, une d’échange et une de corrélation comme suit: xc = x + c (1.15) 10 Chapter 1. Rappels sur la théorie de la Fonctionnelle de la Densité (DFT) Si l’on néglige la contribution de la fonctionnelle de corrélation c, on obtient la contribution de Hartree-Fock, reformulée par Dirac sous la forme suivante : x[ρ] = − 3 4πrs  9π 4 1/3 en Hartree (1.16) avec rs =  3 4πρ1/3 1 a0 , a0 est le rayon de Bohr et 1 Hartree ↔ 27.211 eV . Une amélioration de la LDA qui permit de prendre en compte le spin de l’électron, indispensable pour les systèmes magnétiques, fut ensuite proposée par Barth & Hedin en 1972 [17], ainsi que par Rajagopal & Callaway en 1973 [18, 14]. Il s’agit de l’approximation de la densité locale de spin (LSDA): E LSDA xc = Z ρ(r)xc(ρ↑(r), ρ↓(r))dr (1.17) avec la fonctionnelle d’échange : x(ρ↑, ρ↓) = − 3 2  3 4π 1/3 1 ρ (ρ 4/3 ↑ + ρ 4/3 ↓ ) (1.18) où ρ↑(r) et ρ↓(r) sont les densités locales d’électrons de spin ↑ et de spin ↓. Par la suite, nous utiliserons la fonctionnelle de corrélation LSDA paramétrée par J.P. Perdew et Y. Wang en 1992 [19] à partir du résultat des simulations MonteCarlo de Ceperley et Alder [20]. Dans toute la suite, nous regarderons plus précisément les fonctionnelles qui ont pris historiquement le pas sur la fonctionnelle locale (LDA), à savoir l’approximation du gradient généralisé (GGA pour Generalized Gradient Approximation). La GGA et la fonctionnelle PBE (J.P. Perdew, K. Burke, M. Ernzerhof) Dans l’approximation du gradient généralisé, l’énergie de corrélation dans un volume élémentaire dépend de la densité locale de spin ↑ et de spin ↓ et du gradient de ces densités. Une des premières fonctionnelles GGA élaborées a été la fonctionnelle PW91 [19, 21], qui est une fonctionnelle construite pour répondre au plus grand nombre de configurations possibles. Plus récente, la fonctionnelle PBE essaie d’avoir le moins de paramètres implémentés possible. Avec la fonctionnelle PBE, les auteurs ont essayé de corriger certains problèmes de PW91. En ce sens, ils ont  développé une fonctionnelle plus rapide à l’utilisation et qui ne dépend pas de trop de paramètres expérimentaux ou préalablement calculés. La fonction analytique correspondant à la fonctionnelle PBE est claire et transparente. La non-homogénéité des paramètres entraînait dans PW91 des artefacts dans le potentiel d’échange et de corrélation.1 L’expression analytique correspondant à la fonctionnelle PW91 ne prend pas en compte la bonne densité de probabilité pour les hautes densités. De plus, l’approximation faite dans PW91, qui est de développer le gradient de la densité jusqu’au deuxième ordre pour les densités qui varient lentement ou faiblement, est moins précise dans le cas général que l’approximation LSDA. D’un point de vu formel, l’énergie de corrélation E GGA c , pour la fonctionnelle PBE, est donnée par: E GGA c [ρ↑, ρ↓] ≡ Z d 3 r ρ[ unif c (rs, ξ) + H(rs, ξ, t)] (1.19) H dépend du gradient local de la densité, rs est le rayon local de la sphère de WignerSeitz (en unité a0), ξ est la polarisation en spin et t est le gradient de densité. Ces deux dernières quantités sont définies comme suit: ξ = ρ↑ − ρ↓ ρ (1.20) t = ∇ρ 2φ(ξ)ksρ (1.21) avec φ(ξ) = (1 + ξ) 2/3 + (1 − ξ) 2/3 2 (1.22) ks = q 4kF /πa0 (1.23) où kF est le vecteur d’onde de Fermi du gaz d’électrons uniforme de densité ρ. L’énergie d’échange et de corrélation pour la fonctionnelle PBE est donnée par : E P BE xc = Z d 3 r ρ unif x Fxc(rs, ζ, s) (1.24) avec s = |∇ρ|/2kF ρ. 1Dans l’article [21], les auteurs proposent trois limites asymptotiques de la théorie : la variation lente (∇ρ → 0), la variation rapide (∇ρ → ∞), et le dimensionnement uniforme dans la limite de haute densité (λ → ∞; λ étant un paramètre d’échelle). Cette dernière limite est caractérisée par Ec[ρ↓, ρ↑ → cste, or dans PW91 cette condition n’est pas respectée, d’où la survenue d’instabilités. 12 Chapter 1. Rappels sur la théorie de la Fonctionnelle de la Densité (DFT) Remarques sur la non prise en compte des forces de Van-der-Waals Etant donné son caractère local, les méthodes DFT basées sur des fonctionnelles telles que GGA ou LDA ne prennent pas directement en compte les forces de London2 qui constituent la principale contribution aux forces faibles de Van der Waals. Dans certains cas (par exemple pour les cristaux constitués de feuillets monoatomiques faiblement liés les uns aux autres tels que le graphite ou les cristaux massifs de dichalcogénures de métaux de transition), la non prise en compte de ces forces est un problème qui mène à une sur-estimation de la liaison entre plans atomiques. Nous verrons plus tard que cela ne concerne pas les systèmes étudiés dans cette thèse, constitués de feuillets atomiques isolés liés de façon covalente à leurs substrats métalliques. 1.3 Le code W ien2k W ien2k [22] est un code de calcul de la structure électronique basé sur la DFT, dont nous venons de décrire les principes. Les développeurs de W ien2k ont écrit des procédures d’initialisation et de calcul afin d’avoir une convergence optimale du code en temps et en précision. La première étape est de créer le fichier de structure qui décrit et modélise notre matériau (fichier nommé case.struct). Il faut savoir que la structure sera reproduite périodiquement dans toutes les directions. Pour savoir si nous n’avons pas fait d’erreur à l’écriture du fichier de structure, nous avons la sous-routine nn qui vérifie que les sphères atomiques ne se chevauchent pas. Ensuite, viennent les routines sgroup et symmetry qui déterminent respectivement le groupe d’espace lié à la structure cristalline que l’on vient de former et génère les matrices des opérations de symétrie liées à ce groupe d’espace. La résolution de l’équation de Kohn-Sham passe par un cycle auto-cohérent. Pour cela, nous devons faire l’hypothèse d’une première densité électronique. Cette opération est réalisée par la sous-routine lstart, basée sur un code atomique relativiste écrit à l’origine par Desclaux. lstart génère les densités électroniques des atomes. Pour 2Ces forces ont étés formulées théoriquement par Fritz London en 1937 pour décrire l’interaction entre les dipôles électriques induits dans deux molécules. Elles sont issues du caractère probabiliste de la densité électronique. Le déplacement des électrons étant aléatoire autour d’une molécule, il est fort probable qu’à un moment donné la densité électronique autour de la molécule ne soit plus isotrope. Il en résulte un moment dipolaire induit. Ce moment dipolaire induit peut interagir avec d’autres moments dipolaires induits, ce qui créé une force entre les deux molécules. 

 Le code W ien2k 

 résoudre l’équation de Kohn-Sham, nous avons besoin d’échantillonner la partie irréductible de la zone de Brillouin. Une fois cet échantillonnage réalisé, nous pouvons utiliser un schéma d’intégration de la zone de Brillouin basé sur des tétraèdres [23]. L’échantillonnage de la zone de Brillouin est réalisé par le programme kgen. Une fois les vecteurs d’onde de Bloch déterminés, il ne nous reste plus qu’à générer une première densité de charge (dstart), consistant en la superposition des densités atomiques générées avec lstart. L’étape suivante est le calcul numérique du potentiel de Kohn-Sham à partir de la densité initiale. La commande lapw0 génère le potentiel de Kohn-Sham. Ce potentiel est constitué du potentiel extérieur vext, du potentiel de Hartree vH et du potentiel d’échange et de corrélation vxc. Le potentiel de Hartree est solution de l’équation de Poisson. Il est calculé par développement multipolaire de Fourrier. Le potentiel d’échange et de corrélation est calculé à partir d’un réseau numérique sur grille. Une fois le potentiel total calculé, il nous faut résoudre l’équation de Kohn et Sham (routine lapw1 ). La commande lapw1 nous permet d’avoir les valeurs propres et les vecteurs propres suite à la diagonalisation de la matrice du Hamiltonien, en tenant compte de la matrice de recouvrement. La commande lapw2 calcule le niveau de Fermi EF , ainsi que la contribution à la densité de charge de chaque état occupé et pour chaque vecteur de Bloch. Ensuite, les charges (partielles) dans les sphères sont obtenues par intégration. Les états de coeur sont calculés grâce à la commande lcore. Une fois les densités d’états de coeur, de semi-coeur et de valence obtenues, mixer les additionne afin de produire la nouvelle densité. Si cette densité répond aux exigences de convergence, alors la procédures s’arrête ; sinon, la procédure recommence à partir de lapw0. Le schéma algorithmique de W ien2k est décrit sur la figure 1.1. Le code W ien2k offre la possibilité de prendre en compte les effets liés au couplage spin-orbite. Ces effets ne sont pris en compte qu’à l’intérieur des sphères atomiques, seules régions de l’espace où ces effets sont importants. Dans le code les effets du couplage spin-orbite sur les états de Bloch, sur la structure de bandes et sur l’énergie de l’état fondamental sont pris en compte, via un calcul de perturbations à l’ordre 2. Le code calcule d’abord les fonctions de Bloch et leurs énergies sans ternir compte du couplage spin-orbite, puis l’hamiltonien du couplage spin-orbite est considéré comme une perturbation. On peut alors calculer les corrections aux fonctions d’onde et aux énergies dues à cette perturbation. Ce calcul fait intervenir des éléments de matrice du terme de couplage spin-orbite entre les états occupés et les états inoccupés, jusqu’à une énergie de coupure située bien au delà du niveau de Fermi et choisie par l’utilisateur du code. Nous allons maintenant décrire les fonctions de base utilisées par W ien2k pour résoudre l’équation de Kohn et Sham. 

Les bases de fonctions utilisées par W ien2k 

 La base des Ondes Planes Augmentées (APW)

 Bien que peu utilisée de nos jours, la base APW est une bonne entrée en matière. Elle permet, entre autres, de fixer certaines idées que nous retrouverons pour des bases plus sophistiquées. Pour commencer, quelle que soit la base utilisée, nous modéliserons le cristal avec un potentiel muffin-tin, où chaque atome est situé dans une sphère atomique (S) centrée sur son noyau. Les sphères sont entourées par une région interstitielle (I). Cette modélisation vient du constat que les électrons proches d’un noyau se comportent comme dans l’atome isolé où ils sont décrits avec des harmoniques sphériques, alors que les électrons loin du noyau sont presque des électrons libres décrits par des ondes planes. L’idée derrière le potentiel muffintin est de mettre des limites entre ce que l’on peut décrire avec des harmoniques sphériques et ce que l’on peut décrire avec des ondes planes. Pour chaque vecteur d’onde de Bloch ~k appartenant à la première zone de Brillouin, on définit les fonctions de base APW par: φ ~k K(~r, E) =    √ 1 V e i(~k+K~ )·~r ~r ∈ I P l,m A α,~k+K~ l,m u α l (r 0 , E) Y l m(rˆ0 ) ~r ∈ S (1.25) Avec V le volume de la cellule élémentaire. La position de chaque électron est repérée par rapport au centre de la sphère atomique α grâce au vecteur r~0 = ~r − ~rα. Nous utiliserons les coordonnées sphériques associées au vecteur r~0 . Pour chaque vecteur d’onde de Bloch ~k, les ondes planes augmentées dépendent des vecteurs du réseau réciproque K~ . Le paramètre A α,~k+K~ l,m doit être déterminé et le paramètre E a la dimension d’une énergie. Les harmoniques sphériques sont les fonctions Y l m, et les fonctions u α l sont les solutions de la partie radiale de l’équation de Schrödinger pour un atome libre et pour l’énergie E. Une bonne connexion entre le milieu interstitiel et la sphère atomique induit que les fonctions APW doivent être continues à la surface des sphères. Pour cela, nous faisons en sorte que les deux parties soient égales en valeur mais pas nécessairement en dérivée normale. Autrement dit, nous souhaitons que les deux régions satisfassent une condition de continuité. Pour l’obtenir, on doit d’abord développer les ondes planes dans la région interstitielle sur une base d’harmoniques sphériques. On obtient alors : 1 √ V e i(~k+K~ )·~r = 1 √ V e i(~k+K~ )· ~rα X l,m i l Jl(| ~k + K~ ||~r|) Y l m ∗ ( ~k + K~ ) Y l m(rˆ0 ) (1.26) Pour satisfaire à la condition de continuité, on pose ~r = R~ α le rayon de la sphère atomique α. Une fois ceci fait, on obtient : A α,~k+K~ l,m = i l e i(~k+K~ )· ~rα Jl(| ~k + K~ ||~r|) u α l (R~ α, E) Y l m ∗ ( ~k + K~ ) (1.27) Les fonctions φ ~k K(~r, E) ne sont pas des fonctions propres (fonctions de Bloch des électrons dans le cristal) mais uniquement des fonctions de base. Elles sont par construction, proches des fonctions propres, ce qui leur confèrent une bonne efficacité. En théorie, on devrait avoir une infinité de termes dans la somme sur l et m qui donne les APW φ ~k K. En pratique, il nous faut cependant tronquer la somme jusqu’à une valeur lmax. Ainsi, les harmoniques sphériques peuvent être nulles jusqu’à 2lmax fois sur un cercle porté par la sphère atomique. Si l’on raisonne en termes de noeuds de la fonction APW, il faut que l’on ait égalité entre le nombre de noeuds par unité de longueur d’un cercle situé sur la sphère atomique des deux cotés de celle-ci, ce qui nous donne : 2lmax 2πRα = 2 2π Kmax (1.28) et donc lmax = RαKmax (1.29) L’équation (1.29) fait correspondre, sur la sphère atomique, les noeuds des ondes planes avec ceux des harmoniques sphériques. Notre calcul sera plus précis si l’on augmente la valeur du produit RKmax et donc de lmax. Cependant, nous ne pouvons pas prendre de valeurs de lmax trop grandes par rapport au produit RKmax. Ceci reviendrait à prendre un polynôme de trop haut degré pour modéliser une fonction passant par quelques points. Une telle configuration mènerait à des instabilités numériques au niveau de la sphère atomique. Les ondes planes auraient trop de noeuds pour s’ancrer, ce qui créerait des instabilités liées au partage des noeuds. Par ailleurs, les rayons des atomes pour un système donné ne doivent pas être trop différents, sous peine de ne pas trouver de valeur de RKmax optimale pour tous les atomes. Pour décrire les états propres avec la méthode APW, nous devons faire correspondre, pour un état donné, l’énergie E à la valeur propre de cet état. Pour résoudre ce problème auto-cohérent, nous devons faire une hypothèse sur la valeur de l’énergie E en espérant que celle-ci soit aussi proche que possible de la valeur propre  n ~k . Ceci fait, nous pouvons établir la matrice de l’Hamiltonien et la matrice de recouvrement. En résolvant l’équation séculaire ainsi construite, nous devrions trouver la valeur propre  n ~k pour le paramètre E considéré. Si  n ~k est différent de E, nous devons faire une nouvelle hypothèse sur la valeur de E et recommencer le cycle auto-cohérent jusqu’à ce que  n ~k = E et que le cycle s’arrête alors. Cette procédure est réalisée pour toutes les bandes  (n=1) ~k ,  (n=2) ~k , etc… 1.4.2 La base des Ondes Planes Augmentées Linéarisées Le problème majeur avec la base APW est qu’il faut déterminer le paramètre E, alors inconnu, pour connaître les vecteurs propres. C’est ce problème que la base LAPW (Linearised Augmented Plane Waves) propose d’éviter, en exprimant les fonctions de base sans paramètre E inconnu. Pour cela, il faut linéariser u α l (r 0 , α ~k ) autour d’une énergie fixe E0, par un développement de Taylor, ce qui nous donne : u α l (r 0 , α ~k ) = u α l (r 0 , E0) + (E0 −  n ~k ) ∂uα l (r 0 , E) ∂E      E=E0 + O(E − E0) 2 (1.30) En incluant ce développement au premier ordre dans l’équation 1.25, on a: 18 Chapter 1. Rappels sur la théorie de la Fonctionnelle de la Densité (DFT) φ ~k K(~r) =    √ 1 V e i(~k+K~ )·~r ~r ∈ I P l,m  A α,~k+K~ l,m u α l (r 0 , E0) + B α,~k+K~ l,m u˙ α l (r 0 , E0)  Y l m(rˆ0 ) r~0 ∈ S (1.31) La condition de continuité à la surface de la sphère ne porte plus désormais sur la seule fonction φ ~k K(~r), mais aussi sur sa dérivée normale, ce qui induit un ajustement optimal des ondes planes avec les harmoniques sphériques. L’omission de O(E −E0) 2 est justifiée si l’on prend l’énergie E0 proche du centre de la bande dominante pour un atome donné. Cette opération étant singulière pour chaque atome, nous aurons plusieurs énergies pivot E n 0,l avec n = 0, 1, 2, 3… (la signification du l ne viendra que plus tard, lorsque nous voudrons raffiner cette méthode grâce à des orbitales locales). L’équation (1.31) devient donc : φ ~k K(~r, E) =    √ 1 V e i(~k+K~ )·~r ~r ∈ I P l,m  A α,~k+K~ l,m u α l (r 0 , En 1,l) + B α,~k+K~ l,m u˙ α l (r 0 , En 1,l)  Y l m(rˆ0 ) r~0 ∈ S (1.32) Pour diminuer le temps de calcul, on peut choisir d’augmenter Rmin. Pour conserver un niveau de précision équivalent, il nous faut conserver le même RminKmax, ce qui nous permet de diminuer Kmax. Ces deux actions consistent à diminuer le nombre d’ondes planes nécessaires à décrire l’espace interstitiel plus petit et donc de recourir à plus d’harmoniques sphériques pour décrire le surcroît de volume des sphères atomiques. En diminuant Kmax, nous diminuons la matrice à diagonaliser, ce qui diminue le temps de calcul lié à la diagonalisation.