DOURNAC.ORG
Français  English
Accueil
Astronomie
Sciences
Philosophie
Informatique
Cv

Informatique > Calcul numérique de géodésiques dans un Univers d'Einstein-de Sitter avec la métrique FLRW



  • 1.Introduction
  • 2.Univers d'Einstein–de Sitter
  • 3.Système différentiel
  • 4.Résolution numérique
  • 5.Résultats
  • 6.Résolution avec Matlab
  • 7.Cas particulier - Solution analytique
  • 8.Conclusion

1.Introduction :

Ce travail a été réalisé sous la direction de Jacques Harthong. L'expansion de l'Univers influence la trajectoire des rayons lumineux appelée aussi géodésique. La relativité restreinte permet d'étudier leurs comportement dans un référentiel galiléen, mais cette approche ne suffit pas dans le cas d'une expansion telle qu'elle est décrite dans les modèles du Big-Bang. Il faut donc utiliser les résultats issus de la relativité générale. Nous nous intéressons ici à la détermination de la trajectoire des rayons lumineux dans la métrique de Friedmann-Lemaître-Robertson-Walker (FLRW). Cette métrique s'écrit :

\begin{equation} \text{d}s^{2}=c^{2}\text{d}t^{2}-R(t)^{2}\text{d}l^{2} \label{eq1} \end{equation}

où la partie spatiale $\text{d}l^{2}$ peut prendre la forme suivante :

\begin{equation} \text{d}l^{2}=\frac{\text{d}r^2}{1-kr^2}+r^2(\text{d}\theta^2+\text{sin}^2\theta\,\text{d}\phi^2) \label{eq2} \end{equation}

avec $k$ le paramètre de courbure. On parlera alors d'espace de type sphérique ($k=1$), plat ($k=0$) ou hyperbolique ($k=-1$).

La figure suivante nous permet de voir le résultat recherché, c'est-à-dire le calcul de la géodésique lumineuse émise à l'instant présent $t_{0}$ par la galaxie de gauche et reçue par la nôtre à l'instant $t_{0}$.


Figure 1 : Trajectoire d'un rayon lumineux dans un espace en expansion

2.Cas d'un espace d'Einstein-de Sitter ($k=0$ et $\Lambda=0$) :

On modélise ici les géodésiques dans un espace plat ($k=0$) avec une constante cosmologique nulle ($\Lambda=0$). Dans ce cas, le facteur d'échelle a une forme analytique (avec $H_{0}$ la constante de Hubble au moment présent, c'est-à-dire à $t=t_{0}$) :

\begin{equation} R(t)=R_{0}\,\bigg[\frac{3\,H_{0}}{2}\,t\bigg]^{2/3} \label{eq3} \end{equation}

Nous pouvons établir le système différentiel régissant le comportement des rayons lumineux. Celui-ci sera résolu numériquement par la méthode de Runge-Kutta. En calcul tensoriel, l'équation vérifiée par les géodésiques est la suivante (avec la convention de sommation d'Einstein) :

\begin{equation} \dfrac{\text{d}^2 u^{i}}{\text{d}s^2}+\Gamma^{i}_{jk}\,\dfrac{\text{d}u^{j}}{\text{d}s}\,\dfrac{\text{d}u^{k}}{\text{d}s}=0 \label{eq4} \end{equation}

La variable "$s$" est un paramètre affine (comme l'abscisse curviligne), à ne pas confondre avec l'intervalle d'espace-temps défini au début. Puisque la matrice du tenseur métrique est diagonale ($g_{ij}=0$ si $i\neq j$), les symboles de Christoffel de seconde espèce sont :

\begin{equation} \left.\begin{align} \Gamma^{i}_{ij}&=\Gamma^{i}_{ji}=\dfrac{1}{2}\,\partial_{j}\,\ln|g_{ii}|\\ \Gamma^{i}_{jj}&=-\dfrac{1}{2\,g_{ii}}\,\partial_{i}\,g_{jj}\,\,\, \text{avec}\,\, i \ne j \end{align}\right. \label{eq5} \end{equation}

Tous les autres symboles $\Gamma^{i}_{jk}$ sont nuls.

3.Système différentiel :

Après le calcul de ces symboles, on obtient donc le système différentiel non-linéaire suivant:


\begin{eqnarray} \begin{aligned} \dfrac{\text{d}^{2}ct}{\text{d}s^{2}} & = -\dfrac{2}{3}\,R_{0}^{2}\,\bigg(\dfrac{3\,H_{0}}{2\,c}\bigg)^{4/3}\,(ct)^{1/3}\,\bigg(\bigg(\dfrac{\text{d}r}{\text{d}s}\bigg)^{2}+r^2\,\bigg(\dfrac{\text{d}\,\theta}{\text{d}s}\bigg)^{2}+r^{2}\,\text{sin}^{2}(\theta)\,\bigg(\dfrac{\text{d}\phi}{\text{d}s}\bigg)^{2}\bigg)\\ \dfrac{\text{d}^{2}r}{ds^{2}} & = -\dfrac{4}{3\,ct}\,\dfrac{\text{d}ct}{\text{d}s}\,\dfrac{\text{d}r}{\text{d}s}+r\,\bigg(\bigg(\dfrac{\text{d}\theta}{\text{d}s}\bigg)^{2}+\text{sin}^{2}(\theta)\,\bigg(\dfrac{\text{d}\phi}{\text{d}s}\bigg)^{2}\bigg)\\ \dfrac{\text{d}^{2}\theta}{\text{d}s^{2}} & = -\dfrac{2}{r}\,\dfrac{\text{d}r}{\text{d}s}\,\dfrac{\text{d}\theta}{\text{d}s}-\dfrac{4}{3\,ct}\,\dfrac{\text{d}ct}{\text{d}s}\,\dfrac{\text{d}\theta}{\text{d}s}+\dfrac{\text{sin}(2\,\theta)}{2}\,\bigg(\dfrac{\text{d}\phi}{\text{d}s}\bigg)^{2}\\ \dfrac{\text{d}^{2}\phi}{\text{d}s^{2}} & = -2\,\dfrac{\text{d}\phi}{\text{d}s}\,\bigg(\dfrac{2}{3\,ct}\,\dfrac{\text{d}ct}{\text{d}s}+\text{cotan}(\theta)\,\dfrac{\text{d}\theta}{\text{d}s}+\dfrac{1}{r}\,\dfrac{\text{d}r}{\text{d}s}\bigg) \end{aligned} \label{eq6} \end{eqnarray}

Afin de s'affranchir de la valeur de $R_{0}$, nous prenons dans ce système différentiel la variable réduite $r'=R_{0}r$, ce qui nous donne le nouveau système :


\begin{eqnarray} \begin{aligned} \dfrac{\text{d}^{2}ct}{\text{d}s^{2}} & = -\dfrac{2}{3}\,\bigg(\dfrac{3\,H_{0}}{2\,c}\bigg)^{4/3}\,(ct)^{1/3}\,\bigg(\bigg(\dfrac{\text{d}r'}{\text{d}s}\bigg)^{2}+r'^2\,\bigg(\dfrac{\text{d}\,\theta}{\text{d}s}\bigg)^{2}+r'^{2}\,\text{sin}^{2}(\theta)\,\bigg(\dfrac{\text{d}\phi}{\text{d}s}\bigg)^{2}\bigg)\\ \dfrac{\text{d}^{2}r'}{ds^{2}} & = -\dfrac{4}{3\,ct}\,\dfrac{\text{d}ct}{\text{d}s}\,\dfrac{\text{d}r'}{\text{d}s}+r'\,\bigg(\bigg(\dfrac{\text{d}\theta}{\text{d}s}\bigg)^{2}+\text{sin}^{2}(\theta)\,\bigg(\dfrac{\text{d}\phi}{\text{d}s}\bigg)^{2}\bigg)\\ \dfrac{\text{d}^{2}\theta}{\text{d}s^{2}} & = -\dfrac{2}{r'}\,\dfrac{\text{d}r'}{\text{d}s}\,\dfrac{\text{d}\theta}{\text{d}s}-\dfrac{4}{3\,ct}\,\dfrac{\text{d}ct}{\text{d}s}\,\dfrac{\text{d}\theta}{\text{d}s}+\dfrac{\text{sin}(2\,\theta)}{2}\,\bigg(\dfrac{\text{d}\phi}{\text{d}s}\bigg)^{2}\\ \dfrac{\text{d}^{2}\phi}{\text{d}s^{2}} & = -2\,\dfrac{\text{d}\phi}{\text{d}s}\,\bigg(\dfrac{2}{3\,ct}\,\dfrac{\text{d}ct}{\text{d}s}+\text{cotan}(\theta)\,\dfrac{\text{d}\theta}{\text{d}s}+\dfrac{1}{r'}\,\dfrac{\text{d}r'}{\text{d}s}\bigg) \end{aligned} \label{eq7} \end{eqnarray}

Une fois la résolution numérique effectuée, nous multiplierons la variable $r_{i}=r'_{i}/R_{0}$ par le facteur d'échelle $R(t_{i})$ afin d'obtenir la distance physique $D_{\varphi}$ du rayon lumineux par rapport à notre galaxie. La valeur de cette distance au temps $t=t_{i}$ est égale à : \begin{equation} D_{\varphi,i}=D_{\varphi}(t=t_{i})=R(t_{i})r_{i}=r'_{i}\,\bigg[\frac{3\,H_{0}}{2c}\,ct_{i}\bigg]^{2/3} \label{eq8} \end{equation}

4.Résolution numérique :

Nous pouvons maintenant résoudre numériquement ce système grâce à la méthode de Runge-Kutta à l'ordre 4 avec pas adaptatif. Le projet initial développé durant ma premiére annèe de l'ENSPS a été codé en langage C. Nous nous servirons plus tard des solveurs "ode" de Matlab pour faciliter l'implémentation de cette modélisation.

Voici le source et le script gnuplot associé:

  •    main.c

  •    geodesics.gp

La variable sphérique "$r'$" étant par définition positive, nous nous assurons de mettre la valeur absolue à chacune de ses utilisations dans l'algorithme. Une fois le code compilé (avec "gcc -lm main.c"), nous passons à l'exécution. Le programme requiert un paramètre qui est la distance de laquelle part le rayon lumineux par rapport à notre référentiel qu'est notre galaxie.

Concernant les conditions initiales dans l'algorithme de Runge-Kutta, nous nous apercevons que les dérivées initiales sont liées entre elles par la métrique de FLRW définie en \eqref{eq1}. On peut ainsi écrire :

\begin{equation} \bigg(\dfrac{\text{d}ct}{\text{d}s}\bigg)_{0}^{2}=\bigg[\dfrac{R(t_{0})}{R_{0}}\bigg]^{2}\bigg[\bigg(\dfrac{\text{d}r'}{\text{d}s}\bigg)_{0}^{2}+r_{0}'^{2}\bigg[\bigg(\dfrac{\text{d}\theta}{\text{d}s}\bigg)_{0}^{2}+\text{sin}^{2}(\theta_{0})\bigg(\dfrac{\text{d}\phi}{\text{d}s}\bigg)_{0}^{2}\bigg]\bigg] \label{eq9} \end{equation}

Cette équation apparaît dans le source "main.c" sous la forme :

wp=powl((3*Ho/(2*cv)*w),2./3)*sqrtl(powl(xp,2)+powl(x,2)*(powl(yp,2)+powl(sinl(y),2)*powl(zp,2)));

avec wp=(d(ct)/ds)0, xp=(dr/ds)0, yp=(dθ/ds)0, zp=(dφ/ds)0 et w=ct0, x=r0, y=θ0, z=φ0. On prend pour l'instant le cas simple où xp=-1 (ceci implique dr<0 quand ds>0), yp=0 et zp=0. A noter que wp ((d(ct)/ds)0) sera choisi positif, ce qui signifie que d(ct)>0 quand ds croît. Nous verrons plus loin le cas où yp≠0. Pour ce qui est de t0, r0, θ0, et φ0, nous choisissons : t0=2/(3H0), r0=3000, θ0=pi/2, φ0=pi/2.

5.Résultats :

Les graphiques suivants ont été produits avec une distance initiale de 3000 Mégaparsecs (1 parsec=3.262 années-lumière=3.1013km), soit environ 9.8 Milliard d'années-lumière.
Une fois l'exécution terminée, nous produisons le graphique geodesics_r.ps ("gnuplot geodesics.gp"). En abscisse est représenté le temps exprimé en Mégaparsec (attention à faire la conversion en années) et en ordonnée la distance toujours en Mégaparsec. La ligne continue est la trajectoire de la galaxie qui s'éloigne de la nôtre (suivant la loi de Hubble) et en pointillé la trajectoire du rayon lumineux arrivant vers nous.


Figure 2 : Trajectoire du rayon lumineux pour (dr/ds)0=-1,(dθ/ds)0=0,(dφ/ds)0=0

Nous voyons que le rayon arrive dans notre galaxie au bout de 7030 Mpc, soit 22.918 milliards d'années. On peut remarquer l'incurvation, ici la concavité, de la trajectoire due à la décélération au fil du temps de l'expansion dans le cas d'un univers d'Einstein-de-Sitter : la lumière a de moins en moins de difficultés à lutter contre cette expansion; comparée à une expansion linéaire, elle mettra donc, de manière progressive tout au long de sa trajectoire, de moins en moins de temps pour nous parvenir . Par la suite, la variable "r" se met à augmenter, ce qui traduit l'éloignement du rayon après avoir atteint l'origine de notre référentiel. La trajectoire est modifiée lorsque nous prenons (dθ/ds)0≠0 :


Figure 3 : Trajectoire du rayon lumineux pour (dr/ds)0=-1,(dθ/ds)0=0.0001,(dφ/ds)0=0

On s'aperçoit que le rayon se rapproche de nous dans un premier temps puis s'éloigne. Nous en concluons qu'une légère déviation initiale ((dθ/ds)0=0.0001) ne permet pas à la lumière d'arriver dans notre référentiel, ce qui peut se comprendre de manière intuitive. Ceci est confirmé par la figure 4 où est représenté θ en fonction de ct:


Figure 4 : Valeur de θ en fonction de (ct) pour (dr/ds)0=-1,(dθ/ds)0=0.0001,(dφ/ds)0=0

θ varie de π/2 à 3π/2 mettant en évidence le passage de la lumière au voisinage de notre galaxie. Regardons maintenant le résultat obtenu en prenant (dθ/ds)0=0 et (dφ/ds)0≠0 :


Figure 5 : Trajectoire du rayon lumineux pour (dr/ds)0=-1,(dθ/ds)0=0,(dφ/ds)0=0.0001

La géodésique obtenue est la même que celle de la figure(3). On peut se rendre compte en effet qu'il y a une symétrie entre θ et φ dans le système différentiel quand on prend un de ces deux paramètres constant et l'autre variable et vice-versa. Ceci est validé par la courbe suivante qui est la même que celle montrée précédemment pour θ (cf Fig(4)).


Figure 6 : Valeur de φ en fonction de (ct) pour (dr/ds)0=-1,(dθ/ds)0=0,(dφ/ds)0=0.0001

Ce résultat confirme aussi qu'il n'y a aucune direction privilégiée dans la propagation de la lumière, autrement dit on vérifie l'isotropie du rayonnement dans cette métrique.

6.Résolution avec Matlab :

Portons à présent le code en Matlab qui dispose de solveurs déjà implémentés.

Le source Matlab est disponible ici : light_geodesics.m

Les résultats obtenus prédémment en langage C sont validés par le programme Matlab. Le solveur utilisé "ode45" est similaire à l'algorithme de Runge-Kutta avec pas adaptatif.


Figure 7 : Trajectoire du rayon lumineux pour (dr/ds)0=-1,(dθ/ds)0=0,(dφ/ds)0=0

7.Cas particulier - Solution analytique :

Dans le cas particulier où nous choisissons (dθ/dt)=(dφ/dt)=0, de sorte que les rayons lumineux qui nous intéressent suivent des trajectoires radiales, nous n'avons pas à manipuler des intervalles angulaires et cette approche ne considérera que deux variables: $r$ et $t$. Puisque pour des géodésiques, $\text{d}s=0$, on peut écrire :

\begin{equation} c^{2}\text{d}t^{2}=R(t)^{2}\text{d}r^{2} \label{eq10} \end{equation} \begin{equation} c\text{d}t=-R(t)\text{d}r \label{eq11} \end{equation}

Nous mettons le signe moins pour bien spécifier que le rayon part dans notre direction. Il vient :

\begin{equation} {\large\int}_{t_{0}}^{t}\frac{c\,\text{d}t}{R(t)}=-{\large\int}_{r_{0}}^{r}\text{d}r \label{eq12} \end{equation} \begin{equation} \Rightarrow\,\,r=r_{0}-\dfrac{3}{R_{0}}\bigg[\dfrac{3\,H_{0}}{2\,c}\bigg]^{-2/3}\,\big[(ct)^{1/3}-(ct_{0})^{1/3}\big] \label{eq13} \end{equation}

ainsi en multipliant par le facteur $R(t)$, on obtient:

\begin{equation} D_{\varphi}(t)=R(t)r=\bigg[\frac{3\,H_{0}}{2}\,t\bigg]^{2/3}R_{0}r_{0}-3(ct)^{2/3}\big[(ct)^{1/3}-(ct_{0})^{1/3}\big] \label{eq14} \end{equation}

soit en considérant la distance comobile initiale $D_{\text{comobile}}=R_{0}r_{0}$ :

\begin{equation} D_{\varphi}(t)=R(t)r=\bigg[\frac{3\,H_{0}}{2}\,t\bigg]^{2/3}\,D_{\text{comobile}}-3(ct)^{2/3}\big[(ct)^{1/3}-(ct_{0})^{1/3}\big] \label{eq15} \end{equation}

nous avons désormais la solution analytique $D_{\varphi}(t)$ que nous traçons, en prenant $D_{\text{comobile}}=3000\,\text{Mpc}$ :


Figure 8 : Trajectoire du rayon lumineux pour (dθ/dt)=0,(dφ/dt)=0

Cette trajectoire est identique à la solution numérique obtenue sur la Figure 2. Ici se trouve donc la solution particulière du problème global explicité depuis le début.

8.Conclusion :

Nous avons ainsi pu étudier la propagation de la lumière dans un espace évoluant avec le temps. Il faut noter que la courbe obtenue pour la gédésique uniquement radiale est dans notre cas de nature concave alors que l'on s'attend, dans le modèle de concordance actuel (ΛCDM), à une courbe convexe (comme ceci est illustré sur la Figure 1) étant donné l'accélération de l'expansion actuellement observée. Dans ce modèle standard, il n'y a pas de formule analytique pour le facteur d'échelle $R(t)$, il faudra donc recourir à d'autres méthodes numériques adaptées à ce problème. Enfin, la mise en évidence du comportement non linéaire des rayons remet en question notre perception euclidienne des phénomèmes physiques à une échelle plus modeste que les distances que nous considérons dans cette étude.


ps : contribuez comme moi au projet Cosmology@Home dont le but est d'affiner le modèle décrivant le mieux notre Univers.

    Accueil | Astronomie | Sciences | Philosophie | Informatique | Cv    
- dournac.org © 2003 by fab -

Haut de Page