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 Friedmann-Lemaître-Robertson-Walker (FLRW) metric. This metric can be written :

$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}$.

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}$) :

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})$) :

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 :

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 :

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 :

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 FLRW metric defined in \eqref{eq1}. Then, we can write :

with wp=(d(ct)/ds)_{0}, xp=(dr/ds)_{0}, yp=(dθ/ds)_{0}, zp=(dφ/ds)_{0}
and w=ct_{0}, x=r_{0}, y=θ_{0}, z=φ_{0}.
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)_{0}) 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_{0}, r_{0}, θ_{0}, and φ_{0}, we take : t_{0}=2/(3H_{0}), r_{0}=3000, θ_{0}=pi/2, φ_{0}=pi/2.

The following graphs were produced with an initial distance of 3000 Megaparsecs (1 parsec = 3.262 light-years = 3.10^{13}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.

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 :

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

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.

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