N-corps évolutif pour la modélisation photométrique et dynamique des galaxies de type précoce

Code Syer-Tremaine généralisé : Méthode

Description/Algorithme

Le principe est simple : les particules d’un N-corps voient leur poids évoluer au cours du temps afin d’ajuster les observables de la galaxie étudiée. La méthode a été proposée par Syer & Tremaine (1996) (voir 1.2.3). Ces derniers ont proposé une formule de prescription de changement de poids des particules au fur et à mesure de leur intégration sur leur orbite. Pour présenter cette méthode, les auteurs ont ajusté des modèles de masse analytiques de plusieurs types (sphère de Plummer, Plummer triaxial,. . .) avec différents jeux de conditions initiales (sphère homogène, distribution suivant celle du modèle de masse). Les calculs des forces sont effectués soit analytiquement soit numériquement par transformation de Fourier (FFT). Des tests ont été menés afin de déterminer le type de conditions initiales ainsi que la valeur des différents paramètres d’entrée les plus appropriés. Cette première étape menée par Syer & Tremaine nous a permis de savoir dans quelle mesure il était possible de développer un code de modélisation dynamique partant sur ces bases. L’originalité du travail présent est d’avoir pu généraliser cette méthode à d’autres types d’observables. En plus d’ajuster la luminosité des galaxies observées, on reproduit leur dynamique, comme suggéré par Bissantz et al. (2004). L’utilisation d’observations de galaxies en spectroscopie intégrale de champ nous permet d’avoir des données détaillées pour contraindre le modèle construit. Le code a été conçu de façon à ce qu’il soit possible d’ajuster parallèlement différentes grandeurs linéaires. Le code permet ainsi de modifier la luminosité (masse) de chaque particule. Les vitesses ne peuvent pas être modifiées, mais les particules ayant une vitesse inappropriée vont voir leur poids diminuer. Dans un N-corps classique, les particules subissent l’attraction gravitationnelle de leurs congénères et forment ainsi un système autogravitant. Dans le cadre de notre travail, les particules évoluent dans un potentiel fixe Φ(r) dont la forme est déterminée par la distribution de lumière de la galaxie que l’on veut reproduire. On suppose que Φ admet une solution stable décrite par une fonction de distribution dans l’espace des phases f(r, v). Une fois ce potentiel déterminé, des conditions initiales sont tirées (ie position des particules dans l’espace des phases). Les particules évoluent par la suite dans ce potentiel, la contribution de chaque particule à chaque observable étant calculée à chaque pas de temps. La comparaison des observables du modèle à celles des données permet de réévaluer le poids attribué à chaque particule. L’ajustement du modèle sur les observations s’effectue par un changement de poids régulier de chaque particule au fur et à mesure de leur orbite, jusqu’à la convergence de leur poids vers une valeur finale. Cette pondération est assimilée à la masse de la particule. 

Code Syer-Tremaine généralisé : Méthode

Le code, entièrement modulaire, peut se diviser de la façon suivante : – détermination du potentiel, – tirage des conditions initiales des particules, – intégrateur, – lecture des données que l’on veut ajuster(soit données réelles, soit analytiques) ainsi qu’écriture dans le format du code, – projection du modèle sur ces mêmes observables, – comparaison des observables du modèle et des données, – calcul & application du changement de poids sur chaque particule L’algorithme est schématisé ci-dessous : FIG. 2.1: Algorithme du code. Dans la suite de ce chapitre nous allons détailler chaque module, en commençant par justifier les unités choisies. Il est en effet essentiel qu’elles soient adaptées aux échelles du modèle. 

Unités

Une galaxie a des dimensions de l’ordre de quelques dizaines de kpc et des vitesses de quelques centaines de km s −1 . Les instruments qui nous permettent d’observer en détail ces galaxies sont adaptés à ces échelles. L’instrument SAURON a par exemple un champ de 4100 × 3300 , avec un échantillonnage spatial de 0.9400 . Prenons une galaxie de l’amas Virgo par exemple, à une distance moyenne de 20 Mpc. SAURON permet de couvrir les 3.5 kpc centraux de cette galaxie, chaque pixel couvrant environ 90 pc. Nous souhaitons modéliser de telles galaxies pour lesquelles des observations sont à notre disposition. Pour un code numérique il convient de ne pas manipuler des chiffres ayant de trop grands rapports afin de ne pas perdre de précision par arrondis. Il est donc inadapté de travailler en S.I. Au vu des échelles galactiques, on a choisi les unités de distances en pc et de vitesses en km s −1 , ce qui fixe l’unité de temps : v = d t (2.1) [t] = 3.085677 1013 3.155815 107 = 0.9778 106 ans (2.2) 28 CHAPITRE 2. Code Syer-Tremaine généralisé : Méthode 2.3. Conditions initiales… Si l’on choisit G = 1, il faut que : 6.67210−11 m3 s 2 kg = 1 (2.3) Cela détermine l’unité de masse, qui est de l’ordre de grandeur de quelques masses solaires : 6.672 10−17 3.0856702 1016 m2 pc s 2 kg = 1 232.52 km2 pc s 2 M (2.4) L’unité de masse choisie est alors de [M] = 232.52M. Ce système d’unités est couramment utilisé en modélisation de galaxies : [G] = 1 [r] = pc [v] = km.s−1 [t] = 0.9778 106ans [M] = 232.52M [L] = 232.52L 2.3 Conditions initiales… On appelle conditions initiales la position des particules dans l’espace des phases, ie leur position spatiale et leur vitesse. Leur masse initiale est aussi fixée. Le choix des conditions initiales est crucial. Le code ne permet de changer que les masses des particules et la création/annihilation des particules ne permet pas de radicalement changer la distribution des particules dans l’espace des phases. C’est pourquoi il paraît essentiel de partir de conditions initiales échantillonnant correctement la distribution réelle. Afin d’effectuer des tests sur le code on ajuste des données analytiques dont les paramètres sont connus. Il est alors aisé de tirer des conditions initiales dans ce type de potentiel. Par contre, si l’on veut modéliser une galaxie réelle, on va devoir déterminer le potentiel dans lequel vont évoluer ces particules. Un choix judicieux du potentiel peut nous permettre de tirer des conditions initiales aisément. Lorsque, dans le futur, les particules évolueront dans un potentiel autocohérent, cette méthode de tirage de conditions initiales sera très probablement conservée.

..dans un potentiel analytique.

On utilise des conditions initiales tirées dans des potentiels analytiques simples soit pour effectuer des tests du code, soit comme point de départ pour certaines galaxies réelles dont la distribution ressemble à une distribution analytique connue. Ont été codées la génération des conditions initiales pour une sphère de Plummer, pour une sphère homogène, ainsi que pour une distribution de Miyamoto-Nagai. …de type Plummer Les conditions initiales de Plummer sont inspirées du papier de Aarseth et al. (1974). Les particules sont initialisées à masses égales, soit m = Mtot/N. De la fonction de distribution, on déduit que la masse à l’intérieur d’une sphère de rayon r s’écrit : M(r) = Mr 3 (r 2 + b 2) 3 2 (2.5) 2

Conditions initiales

 Afin de déterminer la position r d’une étoile, on génère simplement un nombre aléatoire X1 1 par lequel on multiplie M la masse totale pour initialiser M(r), on a ainsi : r = b q X − 2 3 1 − 1 (2.6) Le vecteur position de l’étoile doit être déterminé sur une sphère de rayon r, avec une probabilité uniforme. Pour cela on construit un vecteur aléatoire unitaire u, dont la norme est ramenée à r en multipliant par le rapport de leurs normes : r = krk kuk u (2.7) On calcule ensuite le module de vitesse pour la même particule. La vitesse maximale de V à la distance R est la vitesse d’échappement : ve(r) = p 2 |Φ(r)| (2.8) soit ve(r) = s 2 GM p (r 2 + b 2) (2.9) On pose V/V e = q. Aarseth et al. (1974) montre que la distribution de probabilité de q est proportionnelle à : g(q) = q 2 (1 − q 2 ) 7/2 (2.10) Un moyen pratique d’échantillonner q en accord avec cette distribution est fourni par la technique de réjection de von Neumann2 . Les valeurs possibles de q sont entre 0 et 1, et g(q) est toujours inférieur à 0.1. On génère alors deux nombres aléatoires normalisés X4 et X5. Si 0.1X5 < g(X4), on adopte q = X4. Si ce n’est pas le cas, une nouvelle paire de nombres aléatoires est tirée, jusqu’à ce qu’on en trouve une satisfaisant l’inégalité. La distribution des vitesses étant isotrope, on peut utiliser la même méthode pour définir un vecteur vitesse de norme v. On répète cette procédure pour toutes les particules. On a ainsi pour chaque particule sa position dans l’espace des phases, correspondant à une distribution de type Plummer. …de type homogène sphérique Les conditions initiales sont beaucoup plus simples à implémenter : on tire aléatoirement sur les trois dimensions spatiales une valeur comprise dans l’intervalle [−b : +b], où b est le rayon de la sphère homogène. Les particules dont le vecteur position est tel que krk > b sont éliminées. On se retrouve donc avec une distribution de sphère homogène. La densité ρ s’écrit : ρ = 3M 4πb 3 (2.11) Dans la version actuelle du module les vitesses des particules ont été initialisées à zéro (ce qui a pour conséquence un effondrement froid). 1avec xrandom, générateur de nombres aléatoires à distribution de probabilité uniforme ente 0 et 1 2Considérons une densité de probabilité f(x), qui peut être délimitée par une boîte (ou une fonction majorante) entre xmin et xmax et ayant un maximum fmax. L’algorithme suivant génère une série de nombres qui suivent f(x) : – générer un nombre aléatoire x, uniformément distribué entre xmin et xmax – générer un second nombre aléatoire indépendant u uniformément distribué entre 0 et fmax – si u < f(x), accepter x. Sinon, rejeter x et répéter l’opération. 

. Conditions initiales… …de type Miyamoto-Nagai

On utilise la même méthode que pour la distribution de Plummer, à savoir que l’on tire une masse M(R) au hasard, entre 0 et Mtot : M = M(R) Mtot = 1 − s (a + b) 2 R2 + (a + b) 2 (2.12) ce qui permet de remonter à R : R = (a + b) p M(2 − M) 1 − M (2.13) Le tirage de z suit aussi la méthode de réjection de von Neumann, en utilisant une fonction majorante. Les vitesses sont tirées selon une distribution normale de largeur σ telle que : σ 2 = GM(a + √ z 2 + b 2) 2 √ z 2 + b 2 2(a + R2 + (a + 3zb)(a + √ z 2 + b 2) 2) q R2 + (a + √ z 2 + b 2) , (2.14) centrée en 0 pour vr, et en v0 pour vθ tel que : v 2 0 = GMaR2 (a + R2 + (a + 3zb)(a + √ z 2 + b 2) 2) q R2 + (a + √ z 2 + b 2) (2.15) 2.3.2 …dans un potentiel galactique (axisymétrique) Pour modéliser des galaxies dont la distribution n’a pas de forme analytiquement simple, on doit initialiser les particules en reproduisant approximativement les distributions de densité et de vitesse observées. Comme mentionné plus haut on doit partir de conditions initiales proches de la réalité, la distribution des particules dans l’espace des phases ne pouvant être radicalement modifiée. Pour illustrer ce propos nous étudierons dans le chapitre suivant des tentatives d’ajustement d’une distribution de Plummer avec des conditions initiales de sphère homogène et vice-versa. On veut déterminer les conditions initiales et calculer l’accélération des particules dans le potentiel de la galaxie que l’on observe. L’idéal serait de travailler avec un potentiel descriptible analytiquement. C’est pourquoi nous nous sommes basés sur un ajustement MGE (Multi-Gaussian Expansion, Monnet et al. (1992); Emsellem et al. (1994)) qui décompose la distribution de lumière de la galaxie en somme de J gaussiennes bidimensionnelles l’ajustant. Celles-ci sont caractérisées par : – leur intensité maximale Ij , – leur rapport masse/luminosité Γj – ce qui permet de définir le paramètre Pj = Γj × Ij – leur largeur σj – leur aplatissement qj , ou excentricité ej telle que e 2 j = 1 − q 2 j – leur centre (xj , yj ), – leur angle de position PAj . Ce modèle analytique de la distribution de luminosité observée tient compte des effets de convolution de la réponse impulsionnelle de l’instrument (Point Spread Function ou PSF). Connaissant l’inclinaison de la galaxie, il est possible de déprojeter ces gaussiennes. On veut alors tirer des conditions initiales assez proches des observations. Afin de déterminerles positions initiales des particules, on peuple chacune des composantes du modèle selon une distribution gaussienne3 centrée sur (xj , yj ) et de largeurs σj . Ces gaussiennes sont tronquées à un (R, z)max fixé, toutefois sans perte de masse, cette troncature étant basée sur un mécanisme de réjection. 3 ce module est codé en python, afin d’utiliser des outils existants : le module numarray.random_array avec la distribution normal

Intégrateur 

Il n’est cependant pas suffisant de ne déterminer que les positions des particules. Il faut aussi connaître leur vitesse initiale. Une vitesse est choisie aléatoirement selon une distribution gaussienne de moyennes (0, 0, vθ) et de dispersions (σR, σz, σθ).