BmyOGG - GPU Numerical Simulation of Galaxy

DOURNAC.ORG
	Français English

Coding > BmyOGG - Building and modeling your Own GPU Galaxy

1.Introduction :

This part presents my current work on a GPU/OpenCL N-Body code applied to the modeling and numerical simulation of galactic dynamics. Graphical user interface is achieved with Qt framework and Graphics rendering with OpenGL VBO (Vertex Buffer Object). Below on the video, blue particles represent dark matter; white ones are "classic particles", i.e a collection of stars. Initial conditions for coordinates and velocities are based on the Springel-Hernquist model with a massive dark matter halo and a stellar disk. The time step in this simulation equals about 1 million years.

After a short introduction to the main features of these objects, we discuss in section 4 about the different theoretical models and numerical tools needed to perform these simulations on the two main types of galaxies: spiral and elliptical

2.Results :

3.Features :

3.1 Historique :

Les galaxies ne sont connues en tant que mondes à part équivalents à notre Voie Lactée que depuis le début du 20ème siècle. Hubble, en 1924, parvint à résoudre la "nébuleuse" d'Andromède en étoiles, et à prouver ainsi qu'il s'agit d'une autre galaxie et non d'une nébuleuse de gaz comme celle d'Orion qui doit son éclat au rayonnement de l'hydrogène ionisé.

3.2 La Voie Lactée :

Une galaxie est un gigantesque amas d'étoiles, de poussières et de gaz. Notre Voie Lactée contient approximativement 100 milliards d'étoiles. C'est une galaxie spirale, véritable disque mince en rotation, d'un diamètre environ égal à 30 kpc (1 pc = 3.26 a.l) et d'une épaisseur de l'ordre de 1 kpc. La densité moyenne d'étoiles est alors égale à 0.1 étoile par pc³.

Le gaz interstellaire ne représente aujourd'hui qu'une faible fraction de la masse totale des galaxies, de l'ordre de 10% dans les galaxies spirales. Pourtant c'est à partir du gaz que se sont formées et se forment encore actuellement les étoiles. Ce gaz est essentiellement de l'hydrogène, sous forme atomique ou moléculaire ou encore parfois ionisé, mélangé à de l'hélium (25 % en masse), et à d'autres éléments.

3.3 Classification des galaxies :

Les principaux types morphologiques sont représentés par le célèbre "diapason" de Hubble sur la figure ci-dessous; les galaxies sont séparées en trois classes principales: elliptiques (E), spirales (S) et irrégulières (Irr). Les spirales se composent de deux familles, les spirales dites "normales" et les spirales barrées (SB), et se divisent en outre en types Sa, Sb, Sc qui correspondent plutôt à une évolution graduelle qu'à des classes bien distinctes. De gauche à droite sur le diapason de Hubble, le disque inexistant dans les elliptiques, prend de l'importance, de même que la proportion de gaz et d'étoiles jeunes. Les galaxies lenticulaires (SO) sont un intermédiaire entre les spirales et les elliptiques. Environ 2/3 des galaxies sont spirales, l'autre tiers se partageant essentiellement entre 10 % d'elliptiques et 20 % de lenticulaires. Les irrégulières ne constituent que quelques pour cent.

3.4 Galaxies elliptiques :

Les galaxies elliptiques (E) sont des concentrations de forme quasi-sphérique de milliards d'étoiles. Elles nous apparaissent en projection sur le ciel comme des ellipses plus ou moins aplaties. Le rapport des axes (a et b) varie de 1 à 3 seulement, l'ellipticité e définie par e=(a-b)/a (a étant le grand axe) varie donc de 0 à 0.7. Les elliptiques sont classées par un indice qui vaut 10e (de E0 à E7). Les galaxies les plus massives que l'on connaisse sont des galaxies elliptiques, qui sont en grande majorité plus lumineuses que les spirales. Il existe pourtant une classe de galaxies elliptiques naines, en général compagnons de galaxies plus grosses. Les galaxies elliptiques n'ont pas de sous-structures particulières, ce qui entraîne des mouvements d'étoiles aléatoires contrairement aux galaxies spirales. Leur luminosité décroît très régulièrement du centre au bord. Elles ne possèdent pas ou très peu de gaz, mais des étoiles de population relativement vieille, et des amas globulaires, qui ne sont que de très vieux amas d'étoiles. Hubble qualifiait ces galaxies de "précoces", pensant qu'elles étaient le stade antérieur dans l'évolution vers une forme spirale. Aujourd'hui, les astrophysiciens pensent au contraire qu'elles sont le résultat de la fusion de deux galaxies spirales (voir ici une simulation).

3.5 Galaxies spirales :

Les galaxies spirales (S) doivent leur appellation bien sûr aux spectaculaires bras spiraux qui s'enroulent dans un disque très mince ; ces structures sont très aplaties comme le montrent certains spécimens vus par la tranche : leur rapport d'axes peut aller jusqu'à 100 ! Les bras spiraux rejoignent vers le centre le composant le plus brillant, le bulbe, dont la forme est ellipsoïdale et bien moins aplatie que le disque; c'est le bulbe qui est à l'origine de la forme renflée au centre d'une galaxie spirale vue de profil. Dans les galaxies spirales barrées (SB), le noyau est traversé par une barre d'étoiles aux extrémités de laquelle débutent les bras spiraux. Ceux-ci sont en général moins enroulés : on parle de spirale ouverte.

Dans les bras spiraux se forment de nombreuses jeunes étoiles, ce qui leur donne cette couleur bleue caractéristique. A l'inverse, le bulbe contient des étoiles agées (de couleur rouge). Il est admis aujourd'hui que dans la plupart des galaxies, le bulbe contient en son centre un trou noir supermassif de l'ordre de plusieurs millions de masse solaire. Le halo stellaire se compose d'amas globulaires, c'est-à-dire des amas d'étoiles isolées et vieilles. Il représente entre 10 et 30 % de la masse totale du disque.

4.Modélisation :

4.1 Modèle de Springel-Hernquist :

Comme nous l'avons dit dans l'introduction, le modèle de Springel-Hernquist permet de produire des conditions initiales (positions et vitesses) pour 2 types de particules, celles du halo massif de matière noire et celles du disque stellaire. Dans ce modèle, le milieu interstellaire ainsi que le bulbe ne sont pas pris en compte.

Halo de matière noire :

La densité de matière noire s'écrit :

$\displaystyle \rho(r)=\dfrac{M_{DM}}{2\pi}\,\dfrac{a}{r}\,\dfrac{1}{(r+a)^{3}}$ (1)

où $M_{DM}$ est la masse totale de la matière noire. Le paramètre depend de la longueur d'échelle du halo, $r_{s}$ , et de la concentration du halo, $c=r_{200}/r_{s}$ , du profil NFW [Springel et al., 2005] :

$\displaystyle a=r_{s}\sqrt{2[(ln(1+c)-c/(1+c)]}$ (2)

La distribution de masse est décrite par :

$\displaystyle M(r)=\dfrac{M_{DM}\,r^{2}}{(r+a)^{2}}$ (3)

et le potentiel gravitationnel par [Hernquist 1990] :

$\displaystyle \Phi=-\dfrac{GM}{r+a}$ (4)

La vitesse circulaire locale s'écrit :

$\displaystyle v_{c}^{2}(r)=r\,\dfrac{\partial\Phi}{\partial r}=\dfrac{GM\,r}{(r+a)^{2}}$ (5)

La masse totale de la matière noire est calculé en reliant la vitesse virielle à la masse virielle du profil NFW [Springel and White, 1999] :

$\displaystyle M_{200}=\dfrac{v_{200}^{3}}{10GH}$ (6)

et $r_{200}$ est donné par:

$\displaystyle r_{200}=\dfrac{v_{200}}{10H}$ (7)

Nous pouvons donc calculer $r_{s}$ à partir de la concentration et en déduire ainsi le paramètre .
Disque :
Le disque est modélisé avec la densité suivante :

$\displaystyle \rho_{d}(R,z)=\dfrac{\Sigma(R)}{2\,z_{0}}\,sech^{2}(\dfrac{z}{z_{0}})\,\,\,\,\,\,\,\,\,$ avec $\displaystyle \,\,\,\,\,sech(z)=(\cosh(z))^{-1}=\dfrac{2}{e^{z}+e^{-z}}$ (8)

et

$\displaystyle \Sigma(R)=\Sigma_{0}\,e^{-R/R_{d}}=\dfrac{M_{d}}{2\pi\,R_{d}^{2}}\,e^{-R/R_{d}}$ (9)

$M_{d}$ est masse du disque, $R_{d}$ est son échelle de longueur du disque, et $z_{0}$ son épaisseur.

Nous supposons que la vitesse locale circulaire correspond à la vitesse circulaire d'un disque purement exponentiel :

$\displaystyle v_{c}^{2}(R)=R\,\dfrac{\partial\Phi}{\partial R}=4\pi\,G\,\Sigma_... ...\,\bigg(\dfrac{R}{2\,R_{d}}\bigg)\,K_{1}\,\bigg(\dfrac{R}{2\,R_{d}}\bigg)\bigg]$ (10)

où et sont les fonctions de Bessel cylindriques régulières respectivement de première et seconde espèce.

Le potentiel du disque est calculé grâce à une méthode de maillage dont la routine est issue du travail mené par [Hockney and Eastwood, 1981].
Combinaison du Disque et du Halo :
Nous utilisons pour cette combinaison les travaux analytiques de [Springel and White, 1999] et de [Mo et al., 1998].

Dans ce modèle, 4 hypothèses sont faites sur la nature du disque. Tous les disques sont supposés être exponentiels et dynamiquement stables par rapport à la formation de barres. Les 2 autres hypothèses concernent le moment cinétique et la masse du halo et du disque :

$\displaystyle J_{d}=j_{d}\,J_{DM}$ (11)

$\displaystyle M_{d}=m_{d}\,M_{DM}$ (12)

òu $j_{d}$ et $m_{d}$ sont les rapports du disque par rapport aux quantités globales. [Mo et al., 1998] définit aussi le moment cinétique global du halo par rapport au paramètre de spin :

$\displaystyle \lambda=\dfrac{J_{DM}\,\vert E\vert^{1/2}}{G\,M_{DM}^{5/2}}$ (13)

avec, en supposant que les particules sont toutes sur des orbites circulaires :

$\displaystyle E=-\dfrac{G\,M_{DM}^{2}}{2\,r_{200}}\,f_{c}$ (14)

et :

$\displaystyle f_{c}=\dfrac{c}{2}\,\dfrac{1-1/(+c)^{2}-2\,(ln(1+c))/(1+c)}{[c/(1+c)-ln(1+c)]^{2}}$ (15)

En utilisant la définition ci-dessus de l'énergie, du paramètre de spin et du moment cinétique, on peut écrire :

$\displaystyle J_{d}=M_{d}\,R_{d}\,v_{200}\,\int_{0}^{r_{200}/R_{d}}\,e^{-u}\,u^{2}\,\dfrac{v_{c}(R_{d}\,u)}{v_{200}}\,du$ (16)

L'équation suivante définit la longueur déchelle pour le disque :

$\displaystyle R_{d}=\dfrac{1}{\sqrt{2}}\,\bigg(\dfrac{j_{d}}{m_{d}}\bigg)\,\lambda\,r_{200}\,f_{c}^{-1/2}\,f_{R}(\lambda,c,m_{d},j_{d})$ (17)

avec

$\displaystyle f_{R}(\lambda,c,m_{d},j_{d})=2\bigg[\int_{0}^{\infty}\,e^{-u}\,u^{2}\,\dfrac{v_{c}(R_{d}\,u)}{v_{200}}\,du\bigg]^{-1}$ (18)

La vitesse circulaire du système est :

$\displaystyle v_{c}(r)=\sqrt{v_{c,d}(r)^{2}+v_{c,DM}(r)^{2}}$ (19)

Dans le cas d'un halo adiabatiquement comprimé, $v_{c,DM}$ s'écrit :

$\displaystyle v_{c,DM}=\sqrt{G(M_{f}(r)-M_{d}(r))/r}$ (20)

avec, [Mo et al.,1998] :

$\displaystyle M_{f}(r)=M_{d}(r)+M(r)(1-m_{d})$ (21)

où est la masse totale contenue à l'intérieur du rayon r.

Pour trouver la longueur déchelle du disque, nous utilisons dans le code la fonction de fit suivante pour $f_{r}$ [Mo et al.,1998] :

$\displaystyle f_{R}\sim\bigg(\dfrac{(j_{d}/m_{d})\,\lambda}{0.1}\bigg)^{-0.06+2... ...(j_{d}/m_{d})\lambda]}\,(1-3m_{d}+5.2m_{d}^{2})\,(1-0.019c+0.00025c^{2}+0.52/c)$ (22)

Connaissant maintenant la longueur d'échelle, le paramètre du halo et les autres paramètres définis par l'utilisateur, nous pouvons construire un modèle.

4.2 Génération des positions initiales :

Halo de matière noire :

En partant de la définition de la densité volumique $\rho(r)$ , nous pouvons relier la différentielle de probabilité à la différentielle de masse par la normalisation suivante :

$\displaystyle dP=\dfrac{dm}{M_{total,DM}}\,\,\,\,\,\,\,\,\,\,$ avec $\displaystyle \,\,\,\,\,\,\,\,\,\,dm=4\pi\,r^{2}\,\rho(r)\,dr$ (23)

avec $M_{total,DM}$ la masse totale du halo. Prenons l'exemple d'une distribution de densité comme celle d'Hernquist (voir eq(1)).

Figure 1 : Fonction de répartition avec , $v_{200}=160$ km.s $^{-1}$ , kpc

$\includegraphics[height=9.83cm width=13.05cm]{cumulative.eps}$

La fonction de répartition s'écrit :

$\begin{displaymath}\begin{array}{ll} P(R\le r)&=\mathlarger{\int}_{0}^{r}\,\dfra... ...2\,a\,r'}{(r'+a)^{3}}\,dr'=\dfrac{r^{2}}{(r+a)^{2}} \end{array}\end{displaymath}$ (24)

La figure Fig 1 montre cette solution avec , $v_{200}=160$ km.s $^{-1}$ , $r_{200}=v_{200}/10H=225.35$ kpc, $r_{s}=22.53$ kpc et donc $a\backsimeq\,1.73\,r_{s}= 38.97$ kpc.

Nous devons maintenant générer les positions initiales qui suivent la fonction de densité de probabilité (pdf) qui définit $P(R\le r)$ en (24). La fonction de répartition étant inversable, nous pourrions nous servir la méthode de la transformée inverse ; Nous allons cependant utiliser une méthode plus générale et utile pour la suite, celle de la méthode de rejet :

Soit une fonction de densité de probabilités. On suppose qu'il existe une densité de probabilités telle que :

$\displaystyle \exists K > 0, \forall x \in \mathbb{R}\,\, f(x) \le K g(x)$ (25)

Soit alors suivant la loi de densité , et $Y \sim \mathcal U ([0 ; K g(Z)])$ . Alors la variable aléatoire $X = \lbrace Z \vert Y \le f(Z) \rbrace$ suit la loi de densité . L'algorithme a donc la forme suivante :
Algorithme de la méthode de rejet
On suppose f, g et K connus.

On tire Z suivant la loi de g et $Y \sim \mathcal U ([0 ; K g(Z)])$

Si $Y \le f(Z)$ , on garde Z.

Sinon, on tire un nouveau couple jusqu'à ce que $Y \le f(Z)$
Dans la pratique, avec une loi uniforme égale à $g(Z)=\dfrac{1}{b-a}$ , nous prenons $K=(b-a)\,f_{max}$ . On démontre que le pourcentage de cas favorables est donné par :

$\displaystyle P(U\le f(x)/(K\,g(x)))=\dfrac{1}{K}$ (26)

Figure 2 : Histogramme de distribution du halo de matière noire avec 10 $^{6}$ tirages

$\includegraphics[height=9.83cm width=13.05cm]{hist_pos.eps}$

En remarquant pour ce profil que $f_{max}=\dfrac{8}{27a}$ pour , il suffit de prendre pour :

$\displaystyle K=r_{200}\,f_{max}=225\,\times\,0.0076=1.71$ (27)

La figure Fig 2 représente la distribution générée (sur 10 $^{6}$ tirages aléatoires). Le nombre de cas favorables est biaisé par la valeur de la fonction de répartition en $r_{200}$ ( $P(R\le r_{200})=0.72$ ). Dans le cas théorique, $P(U\le f(x)/(K\,g(x)))=\dfrac{1}{K}$ ; nous pouvons alors déduire qu'en moyenne, nous obtiendrons $\dfrac{0.72}{K}=\dfrac{0.72}{1.71}=0.42$ de cas favorables, soit 42 %.

Il reste désormais à produire les 2 autres v.a $\theta$ et $\phi$ . Nous les distribuons uniformément respectivement sur [ $0,\pi$ ] et [ $0,2\pi$ ] ; ainsi nous déduisons les coordonnées initiales des particules de matière noire et nous les injectons dans le code.
Disque :
Nous suivons le même raisonnement que précdemment pour les positions initiales des particules du disque. La pdf dépend ici des variables et . On peut écrire :

$\displaystyle dP=\dfrac{dm}{M_{total,d}}\,\,\,\,\,\,\,\,\,\,$ avec $\displaystyle \,\,\,\,\,\,\,\,\,\,dm=2\pi\,r\,dr\,dz\,\rho(r,z)$ (28)

où $\rho(r,z)$ est égal à (8):

$\displaystyle \rho_{d}(r,z)=\dfrac{M_{tot,d}}{\pi\,r_{d}^{2}}\,e^{-r/r_{d}}\,\dfrac{1}{4\,z_{0}}\,sech^{2}(\dfrac{z}{z_{0}})\,\,\,\,\,\,\,\,\,$ avec $\displaystyle \,\,\,\,\,sech(z)=(\cosh(z))^{-1}=\dfrac{2}{e^{z}+e^{-z}}$ (29)

Étant donné que les v.a et sont indépendantes, nous pouvons créer les deux distributions indépendamment l'une de l'autre. On peut reécrire (28) sous la forme d'un produit de deux pdf :

$\displaystyle P(R\le r \cap Z\le z)=\mathlarger{\int}_{0}^{r}\,\dfrac{2\,r}{r_{... ...rger{\int}_{0}^{r}\,pdf_{r}(r)\,dr\,\,\mathlarger{\int}_{0}^{z}\,pdf_{z}(z)\,dz$ (30)

Figure 3 : Histogramme de distribution pour la partie avec 10 $^{6}$ tirages ( $r_{d}=3$ kpc et $r_{disk}=20$ kpc)

$\includegraphics[height=9.83cm width=13.05cm]{hist_disk_posr.eps}$

Pour la pdf dépendant de , nous utlisons le même principe que précédemment, c'est-à-dire la méthode de rejet ; pour cela, nous choisissons pour , avec $r_{d}=3$ kpc et $r_{disk}=20$ kpc : $K=(b-a)\,f_{max}=r_{disk}\,f(r_{d})=4.9051$ . La proportion de cas favorables est donc :

$\displaystyle \dfrac{F(r_{disk})}{K}=\bigg[(b-a)\,f_{max}\bigg]^{-1}\,F(r_{disk... ...{disk})=\bigg[\dfrac{2\,r_{disk}}{r_{d}}\,e^{-1}\bigg]^{-1}\,F(r_{disk})=0.4037$ (31)

Sur la figure 3 est représenté l'histogramme de la distribution produite. Faisons de même pour la pdf dépendante de . En prenant $z_{0}=1$ kpc, on a : $K=z_{0}\,g_{max}=z_{0}\,g(0)=0.25$ , ce qui donne une proportion de succès : $\dfrac{G(z_{0})}{K}=\dfrac{G(z_{0})}{z_{0}\,g(0)}=\dfrac{0.1904}{0.25}=0.7616$ . Nous obtenons la figure 4 :

Figure 4 : Histogramme de distribution pour la partie avec 10 $^{6}$ tirages

$\includegraphics[height=9.83cm width=13.05cm]{hist_disk_posz.eps}$

Figure 5 : Histogramme de la distribution du couple avec 10 $^{6}$ tirages

$\includegraphics[height=9.83cm width=13.05cm]{hist_pos_disk_all_1.eps}$

Combinons le produit des deux pdf qui donne une distribution du couple . Nous obtenons, avec $z_{0}=0.3$ , les deux figure 5 et 6.

Figure 6 : Histogramme de la distribution du couple avec 10 $^{6}$ tirages

$\includegraphics[height=9.83cm width=13.05cm]{hist_pos_disk_all_2.eps}$

ps : join like me the Cosmology@Home project whose aim is to refine the model that best describes our Universe