Computation of light geodesics in Einstein-de Sitter Universe with Robertson-Walker (RW) metric

DOURNAC.ORG
	Français English

Coding > Computation of light geodesics in Einstein-de Sitter Universe with Robertson-Walker (RW) metric

1.Introduction
2.Case of Einstein–de Sitter Universe
3.Differential system
4.Numerical solving
5.Results
6.Computation with Matlab
7.Special case - Analytical solution
8.Conclusion

1.Introduction :

This project has been supervised by Jacques Harthong. The expansion of the Universe affects the trajectory of light rays which is also called a geodesic. Relativity allows us to study their behavior in a Galilean reference frame, but this approach is not sufficient in the case of expansion as described in the models of the Big Bang. So one has to use the results of general relativity. We are interested here in determining the path of light rays in the Robertson-Walker (RW) metric. This metric can be written :

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

where spatial term $\text{d}l^{2}$ can be expressed as :

\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}

$k$ is a constant representing the curvature of the space. This is known as spherical ($k = 1$), flat ($k = 0$) or hyperbolic ($k = -1$) space.

The following figure shows the expected result, i.e the computation of the light geodesic emitted at present time $t_{0}$ by the galaxy on the left and received by ours at the moment $t_{0}$.

**Figure 1 : Trajectory of a light ray in an expanding Universe**

2.Case of Einstein–de Sitter Universe ($k=0$ and $\Lambda=0$) :

We model the geodesics in a flat space ($k = 0$) with a zero cosmological constant ($\Lambda = 0$). In this case, the scale factor has an analytical form at present time $t=t_{0}$ (caution, $H_{0}$ is the current Hubble constant, i.e at $t=t_{0}$) :

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

We can establish the differential system governing the behavior of light rays. It will be solved numerically by the Runge-Kutta method. In tensor calculus, the equation satisfied by geodesics is as follows (with the Einstein summation convention and the system of curvilinear coordinates $(u^{i})$) :

\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}

The variable "$s$" is a scalar parameter of motion (like curvilinear abscissa), not to be confused with the space-time interval defined at the beginning. Since the matrix of the metric tensor is diagonal ($g_{ij} = 0$ if $i \neq j$), the Christoffel symbols of second kind are :

\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{with}\,\, i \ne j \end{align}\right. \label{eq5} \end{equation}

All other symbols $\Gamma^{i}_{jk}$ are zero.

3.Differential system :

After calculating these symbols, we get the following nonlinear differential system :

\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}

In order to get rid of the value of scale factor $R_{0}$, one takes the reduced variable $r'=R_{0}r$, that gives the new differential system :

\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}

Once numerical solving is performed, we multiply the discrete variable $r_{i}=r'_{i}/R_{0}$ by scale factor $R(t_{i})$ to get the physical distance $D_{\varphi}$ of light ray relatively to our galaxy. This distance at time $t=t_{i}$ is expressed as :

\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.Numerical solving :

We can now solve this system numerically by the Runge-Kutta order 4 with adaptive step. The initial project developed during my first year at ENSPS was coded in C language. We will use later Matlab "ode" solvers to make easier the implementation of this model.

Here are the source and the gnuplot script for plotting :

main.c

geodesics.gp

Since spherical variable "$r$" is positive by definition, we make sure to set absolute value for each of its uses in the algorithm. Once the code is compiled (with "gcc -lm main.c"), we proceed to execution. The program requires a parameter which is the distance from which the light ray goes away relatively to our galaxy, i.e our inertial frame.

About initial conditions in Runge-Kutta method, we find that initial derivates are related by the RW metric defined in \eqref{eq1}. Then, we can write :

\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}

This equation appears in the source "main.c" as :

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)));

with wp=(d(ct)/ds)₀, xp=(dr/ds)₀, yp=(dθ/ds)₀, zp=(dφ/ds)₀ and w=ct₀, x=r₀, y=θ₀, z=φ₀. Firstly, we take the simple case where xp=-1 (this implies dr<0 when ds>0), yp=0 and zp=0. Notice that wp ((d(ct)/ds)₀) will be chosen as positive, which means that d(ct)>0 when ds grows. We will see later the case with yp≠0. Concerning t₀, r₀, θ₀, and φ₀, we take : t₀=2/(3H₀), r₀=3000, θ₀=pi/2, φ₀=pi/2.

5.Results :

The following graphs were produced with an initial distance of 3000 Megaparsecs (1 parsec = 3.262 light-years = 3.10¹³km), or about 9.8 Billion light-years.
Once the execution is complete, we produce the graph geodesics_r.ps ("gnuplot geodesics.gp"). On x axis is represented local time in Megaparsec (pay attention to the conversion in years) and on y axis distance still in Megaparsec. The solid line is the trajectory of the galaxy moving away from us (from Hubble's law) and dotted the path of the light beam coming towards us.

**Figure 2 : Light ray trajectory with (dr/ds)₀=-1,(dθ/ds)₀=0,(dφ/ds)₀=0**

We see that light ray reach us at Local time equal to 7030 Mpc, i.e 22.918 Billions light-years. We can notice the curvature, here concavity, of the trajectory caused by the deceleration of expansion over time in Einstein-de-Sitter universe : light has less and less difficulty in countering this expansion; compared to a linear expansion, light will take, all over its trajectory, a time increasingly short to come in our galaxy. Subsequently, variable "r" begins to increase, meaning that the light ray, after coming, is going away from our galaxy. Trajectory is modified when we take (dθ/ds)₀≠0 :

**Figure 3 : Light ray trajectory with (dr/ds)₀=-1,(dθ/ds)₀=0.0001,(dφ/ds)₀=0**

One can see that light ray firstly approaches us and then goes away. We conclude that a slight initial deviation ((dθ/ds)₀=0.0001) does not allow light ray to come in our inertial reference, which can be understood intuitively. This is confirmed by Figure 4, which represents θ as a function of ct :

**Figure 4 : θ value as a function of (ct) with (dr/ds)₀=-1,(dθ/ds)₀=0.0001,(dφ/ds)₀=0**

θ varies from π/2 to 3π/2 highlighting the path of light near to our galaxy. Now look at the result obtained by taking (dθ/ds)₀=0 and (dφ/ds)₀≠0 :

**Figure 5 : Light ray trajectory with (dr/ds)₀=-1,(dθ/ds)₀=0,(dφ/ds)₀=0.0001**

This geodesic is the same as that of the figure (3). Indeed, one can see that there is a symmetry between θ and φ in the differential system. You can take one of these two parameters constant and the other variable, and vice versa. This is validated by the following plot wich is the same as that shown previously for θ (see Fig (4)).

**Figure 6 : φ value as a function of (ct) with (dr/ds)₀=-1,(dθ/ds)₀=0,(dφ/ds)₀=0.0001**

This result also confirms that there is no preferred direction in the propagation of light, i.e it satisfies the isotropy of the radiation in this metric.

6.Computation with Matlab :

Let's use Matlab which has already got implemented solvers.

Matlab source is dwonloadable here :

light_geodesics.m

Previous results in C are validated by the Matlab program. The solver used "ode45" is similar to the Runge-Kutta algorithm with adaptive step.

**Figure 7 : Light ray trajectory with (dr/ds)₀=-1,(dθ/ds)₀=0,(dφ/ds)₀=0**

7.Special case - Analytical solution :

In the special case where we choose (dθ/dt)=(dφ/dt)=0, so that the light rays that interest us follow radial trajectories, we have not to manipulate angular intervals and this approach only considers two variables: $r$ and $t$. As for geodesics, $\text{d}s = 0$, we can write :

\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}

We put the minus sign to clearly specify the light ray goes away in our direction. It comes :

\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}

then by multiplying by scale factor $R(t)$, one gets :

\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}

so considering initial comoving distance $D_{\text{comoving}}=R_{0}r_{0}$ :

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

finally, we have the analytical expression of $D_{\varphi}(t)$ that we plot, with $D_{\text{comoving}}=3000\,\text{Mpc}$ :

**Figure 8 : Light ray trajectory with (dθ/dt)=0,(dφ/dt)=0**

This solution is the same as the numerical solution shown on Figure 2. So here is a particular solution of the overall problem explained from the beginning.

8.Conclusion :

We were able to study the propagation of light rays in a space evolving over time. One has to notice that radial trajectory of geodesics obtained has a concave nature whereas, into concordance ΛCDM model, one expects to get a convex curve (as illustrated on Figure 1), given the acceleration of expansion currently observed. In this standard model, there is no analytical expression for scale factor $R(t)$, so we have to use other numerical methods more adapted to this issue. Finally, highlighting the nonlinear behavior of these beams challenges our perception of Euclidean phenomena at physical scale smaller than the distances we have considered in this study.

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