Modélisation le code Particle-in-Cell (PIC)

Le domaine de la simulation des plasmas est généralement compartimenté en deux approches. Lorsque les interactions coulombiennes dominent devant les effets cinétiques (si la densité d’énergie potentielle électrostatique est supérieure à la densité d’énergie cinétique), le plasma est dit couplé, ou dense, et présente des propriétés proches d’un fluide. Pour de très petits échantillons, il est possible de traiter les N2 interactions nécessaires pour résoudre le problème à N corps ; c’est le cas dans les codes de dynamique moléculaire. Lorsque la taille de l’échantillon augmente, cette méthode devient trop coûteuse, et on résout alors les équations fluides du plasma. Dans le cas contraire, si les effets cinétiques sont prépondérants devant les interactions coulombiennes, le plasma est dit cinétique. Un grand nombre de données est alors nécessaire pour résoudre correctement le problème : on résout en fait la fonction de distribution de chaque espèce de particules qui est définie dans un espace à 6 dimensions (positions et impulsions). Les codes qui traitent de ce genre de plasmas sont les codes Vlasov et les codes particulaires. Les fonctions de distribution des différentes espèces sont discrétisées sur un maillage dans les codes Vlasov, et sur leurs caractéristiques dans les codes particulaires. Dans le cas des plasmas considérés en accélération par sillage laser, les codes utilisés appartiennent à cette dernière catégorie, et sont appelés codes PIC (pour Particle-In-Cell). Dans cette partie, nous allons tout d’abord présenter les caractéristiques générales des deux codes PIC qui ont été utilisés dans ce manuscrit, et que sont CALDER [Lefebvre et al., 2003] et CALDER-Circ [Lifschitz et al., 2009]. Ensuite, nous décrirons et validerons la méthode de calcul du rayonnement. III.1 Fonctionnement d’un code PIC III.1.1 Principe général de la méthode PIC Prenons des ordres de grandeur typiques d’une simulation d’accélération par sillage laser. Afin de pouvoir prendre en compte la formation du sillage ainsi que la distribution transverse du laser, et en supposant une densité de 5 × 1018 cm−3 , le plasma qui doit être simulé peut être assimilé à une boîte carrée d’environ 100 µm de côté, ce qui représente N = 5 × 1012 particules. Le nombre de particules n’est pas seulement trop grand pour pouvoir considérer les N2 interactions entre particules ; il l’est déjà pour simplement traiter les N particules individuellement. Ainsi, plutôt que de considérer les Ni particules de l’espèce i, on considère la fonction de distribution fi pour cette espèce. fi(x, p, t)δV δp3 représente alors le nombre de particules de l’espèce i présentes dans un petit volume δV centré en x et dans un petit volume de l’espace des impulsions δp3 centré en p à un temps t. Dans l’approche cinétique et en négligeant les collisions entre particules, l’évolution de fi est alors donnée par l’équation de Vlasov : ∂fi ∂t + vi∇fi + qi(E + vi × B) ∂fi ∂pi = 0, (93) avec qi la charge de l’espèce i. L’équation de Vlasov décrit l’évolution de la répartition des particules soumises à la force de Lorentz, qui dépend des champs électromagnétiques présents. Afin de pouvoir totalement résoudre le système (particules+champs), elle doit donc être couplée avec les équations de Maxwell (3)-(6). Les équations de Maxwell permettent alors de décrire l’évolution des champs à partir des termes sources que sont la densité de courant J et de charge ρ. Ces densités de charge totale et de courant total sont obtenues en sommant les deux premiers moments des distributions fi pour chaque espèce i : ρ = X i ρi = X i qi Z fi(x, p, t)dp, (94) J = X i Ji = X i qi Z vfi(x, p, t)dp. (95) Toutes ces équations forment le système de Maxwell-Vlasov et suffisent à décrire complètement l’évolution du système (particules+champs). Cependant, la résolution numérique de l’équation de Vlasov sur un maillage des coordonnées x et des impulsions p serait trop coûteuse et limitée à 2 × 2 dimensions. Au lieu de cela, dans l’approche PIC, la résolution de Vlasov se fait en discrétisant la fonction de distribution fi pour chaque espèce : fi(x, p, t) = X Nm j Wjδ(x − xj (t))δ(p − pj (t)). (96) C’est-à-dire que la fonction de distribution est maintenant définie comme une somme finie de diracs de poids Wj et centrés aux coordonnées (xj (t), pj (t)) au temps t, que l’on appelle macro-particules. Une macro-particule regroupe ainsi plusieurs particules ayant des coordonnées très proches dans l’espace des phases 6D et pouvant être approximées par (xj (t), pj (t)). Elle possède un poids Wj qui détermine le nombre d’électrons réels qu’elle représente. On peut alors résoudre l’équation de Vlasov par la méthode des caractéristiques en remarquant que fi se conserve le long des trajectoires (xj (t), pj (t)). En conséquence de l’application de cette méthode, la macro-particule suit une trajectoire (xj (t), pj (t)) identique à celle qu’aurait suivie une particule réelle initialement placée à la même position. Ceci permet d’identifier les trajectoires des macro-particules à des trajectoires physiques. Les équations de Maxwell sont quant-à elles en général résolues à l’aide d’un schéma aux différences finies, sur un maillage discrétisant l’espace. Ceci implique 53 III.1 Fonctionnement d’un code PIC d’interpoler successivement les champs connus aux points de ce maillage vers les positions des macro-particules xj (t), puis de projeter les charges et les courants des macro-particules connues aux positions xj (t) sur les points du maillage. Au final, la résolution numérique du système impose de répéter les opérations suivantes à chaque itération temporelle (cf figure 13) : – interpolation des champs E(t) et B(t) du maillage vers les positions des macro-particules xj (t). – résolution des équations du mouvement des macro-particules pour calculer les nouvelles positions xj (t + ∆t). – projection des densités de charge ρj (t+ ∆t) et de courant Jj (t+ ∆t) vers les nœuds du maillage – résolution des équations de Maxwell pour le calcul des champs E(t + ∆t) et B(t + ∆t). La simplicité de ce schéma rend la méthode PIC assez intuitive à implémenter, ce qui, avec la variété des problèmes physiques qu’elle permet de représenter, explique son succès en physique des plasmas à partir des années 60 [Dawson, 1983, Hockney & Eastwood, 1988, Birdsall & Langdon, 2004]. Après avoir vu le fonctionnement global de la méthode PIC, nous allons dans la suite développer successivement les différentes étapes de la boucle PIC pour introduire différents paramètres numériques importants.

 La boucle PIC

Nous allons maintenant voir le fonctionnement plus précis de la boucle PIC. La résolution des équations de Maxwell impliquant un maillage spatial, on introduit les pas spatiaux ∆x, ∆y et ∆z dans les 3 directions de l’espace, ainsi que le pas de temps ∆t. De plus, on introduit la notation En i,j,k pour les champs, qui désigne le champ E au temps n∆t et au point de maillage (i∆x, j∆y, k∆z). De même, on introduit la notation (x n m, p n m) pour la macro-particule m (où on a abandonné la notation concernant le numéro de l’espèce). La résolution des équations de Maxwell se fait en général en utilisant un schéma aux différences finies de type saute-mouton appelé schéma de Yee [Yee et al., 1966]. Cette méthode de résolution implique de connaître E et B de manière décalée dans le temps d’un demi pas de temps ∆t/2. De plus, les différentes composantes des champs E et B ne sont pas connues en des points spatiaux identiques, mais en des points pouvant être décalés d’une demi maille selon le champ et la direction considérée, comme indiqué en figure 14. De la même manière, les positions des particules et leur impulsions sont connues à des temps séparés de ∆t/2. Ces décalages spatiaux et temporels facilitent l’obtention de schémas centrés d’ordre 2.

Interpolation des champs 

Les équations de Maxwell étant résolues sur un maillage spatial, les champs E et B sont uniquement connus en certains points de ce maillage. Avec les notations définies précédemment, on a un accès direct dans le code aux valeurs calculées aux points du maillage : E n xi+1/2,j,k , E n yi,j+1/2,k , E n zi,j,k+1/2 , B n+1/2 xi,j+1/2,k+1/2 , B n+1/2 yi+1/2,j,k+1/2 et B n+1/2 zi+1/2,j+1/2,k . 

CALDER-Circ et autres codes 

La partie précédente a présenté le fonctionnement général d’un code PIC cartésien tel que CALDER. Afin d’accélérer les calculs, ce code a été massivement parallélisé à l’aide de la bibliothèque MPI, ce qui permet d’utiliser les ressources des super-calculateurs : le code est lancé sur de nombreux processeurs qui se partagent le calcul en découpant le maillage spatial en sous-domaines. Pour les simulations habituelles de sillage laser, on ne simule pas toute l’étendue du plasma à chaque pas de temps, mais on utilise une fenêtre glissante qui suit l’impulsion laser pour minimiser la taille du maillage nécessaire. Cependant, même avec ces techniques et avec les moyens numériques croissants qui sont à disposition, les études complètes en 3 dimensions restent très coûteuses, et restent par exemple inaccessibles pour réaliser des simulations lorsque la densité du plasma est très faible, ce qui implique des accélérations sur de plus grandes distances, ou lorsque la durée du laser considéré devient grande (sachant qu’il faut toujours résoudre la longueur d’onde laser dans la direction de propagation). Pour ces raisons, l’amélioration des codes PIC passe non seulement par l’implémentation d’effets nouveaux (collisions, ionisation, effets d’électrodynamique quantique…), mais aussi par le développement de nouveaux codes réduits permettant d’accélérer le calcul en profitant d’approximations pouvant être réalisées dans certains cas. Nous présentons dans la suite les autres codes qui ont été utilisés au cours de cette thèse. 

Code PIC en géométrie cylindrique : CALDER-Circ 

L’importance de considérer une géométrie 3D pour la modélisation de l’accélération par sillage laser, ainsi que le coût élevé des simulations 3D ont motivé l’élaboration de codes réduits. Pour le sillage laser, il est possible de tirer partie de la symétrie de révolution de l’enveloppe du laser ainsi que des champs plasma présents dans son sillage pour se “rapprocher” d’une géométrie cylindrique et réduire le coût des calculs. Cette approche a été utilisée pour la première fois par [Lifschitz et al., 2009] dans le code CALDER-Circ.