Azimuthal geodesics in closed FLRW cosmological models

In the following, we are interested in the geodesic motion in a Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime with line element

$$\mathrm{d}s^{2}=-\mathrm{d}t^{2}+a^{2}(t)\left[\frac{\mathrm{d}r^{2}}{1-kr^{2}}+r^{2}\mathrm{d}\theta^{2}+r^{2}\sin^{2}\negmedspace\theta\mathrm{d}\varphi^{2}\right],$$

(2.1)

where $k$ can take values $\{-1,0,1\}$.

A constant time slice of this spacetime corresponds to a maximally symmetric three-dimensional space with constant positive ($k=1$), zero ($k=0$) or negative ($k=-1$) curvature.

When $k=1$, the so-called scale factor $a(t)$ evolves according to the Friedmann equations

	$\displaystyle\left(\frac{\dot{a}}{a}\right)^{2}$	$\displaystyle=\frac{8\pi}{3}\left(\rho_{1}+\rho_{2}\right)-\frac{1}{a^{2}}\,,$		(2.2)
	$\displaystyle\frac{\ddot{a}}{a}$	$\displaystyle=-\frac{4\pi}{3}\left[\left(\rho_{1}+3p_{1}\right)+\left(\rho_{2}+3p_{2}\right)\right],$		(2.3)

where the dot denotes the derivative with respect to cosmic time $t$. The evolution of $a(t)$ is determined by the energy-momentum content and we assume the universe to be filled with two perfect, barotropic fluids with equation of state $p_{i}(t)=w_{i}\rho_{i}(t)$, for $i\in\{1,2\}$, where $\rho_{i}(t)$ is the energy density, $p_{i}(t)$ is the pressure, and $w_{i}$ is the equation of state parameter. Equations (2.2) and (2.3) imply the energy-momentum conservation,

$$\dot{\rho}_{i}+3\frac{\dot{a}}{a}(\rho_{i}+p_{i})=0\,,\quad\text{for}\;i\in\{1,2\}\,,$$

(2.4)

which is satisfied by each fluid component, and yields

$$\rho_{i}=\rho_{i}^{(0)}a^{-3(1+w_{i})}\,,\quad\text{for}\;i\in\{1,2\}\,,$$

(2.5)

where $\rho_{i}^{(0)}$, for $i\in\{1,2\}$, is a constant of integration. We now define $\beta_{i}\coloneqq 8\pi\rho_{i}^{(0)}/3$, and $\gamma_{i}\coloneqq 1+3w_{i}$, for $i\in\{1,2\}$, and re-write Eq. (2.2) as

$$\left(\dot{a}\right)^{2}=\beta_{1}a^{-\gamma_{1}}+\beta_{2}a^{-\gamma_{2}}-1\,.$$

(2.6)

The geodesic equation is given by

$$\frac{{\rm d}^{2}{x}^{\mu}}{{\rm d}\lambda^{2}}+\Gamma^{\mu}_{\alpha\beta}\frac{{\rm d}{x}^{\alpha}}{{\rm d}\lambda}\frac{{\rm d}{x}^{\beta}}{{\rm d}\lambda}=0\,,$$

(2.7)

and can be derived by considering the Lagrangian $\mathcal{L}$

$$\mathcal{L}=g_{\mu\nu}\frac{\mathrm{d}x^{\mu}}{\mathrm{d}\lambda}\frac{\mathrm{d}x^{\nu}}{\mathrm{d}\lambda}\,,$$

(2.8)

with affine parameter $\lambda$. Note that $\mathcal{L}=0$ for null geodesics and $\mathcal{L}=-1$ for timelike geodesics, respectively. For the metric (2.1), we have the Lagrangian

$$\mathcal{L}=-{t^{\prime}}^{2}+a^{2}(t)\left[\frac{{r^{\prime}}^{2}}{1-kr^{2}}+r^{2}{\theta^{\prime}}^{2}+r^{2}\sin^{2}\negmedspace\theta{\varphi^{\prime}}^{2}\right],$$

(2.9)

with the prime denoting derivatives with respect to $\lambda$. The four components of the geodesic equation (2.7) are

	$\displaystyle{t^{\prime\prime}}+a\dot{a}\left(\frac{{r^{\prime}}^{2}}{1-kr^{2}}+r^{2}{\theta^{\prime}}^{2}+r^{2}\sin^{2}\theta{\varphi^{\prime}}^{2}\right)$	$\displaystyle=0\,,$		(2.10)
	$\displaystyle{r^{\prime\prime}}+2\frac{\dot{a}}{a}{t^{\prime}}{r^{\prime}}+\frac{kr}{1-kr^{2}}{r^{\prime}}^{2}-r\left(1-kr^{2}\right)\left({\theta^{\prime}}^{2}+\sin^{2}\theta{\varphi^{\prime}}^{2}\right)$	$\displaystyle=0\,,$		(2.11)
	$\displaystyle{\theta^{\prime\prime}}+2\frac{\dot{a}}{a}{t^{\prime}}{\theta^{\prime}}+\frac{2}{r}{r^{\prime}}{\theta^{\prime}}-\cos\theta\sin\theta{\varphi^{\prime}}^{2}$	$\displaystyle=0\,,$		(2.12)
	$\displaystyle{\varphi^{\prime\prime}}+2\frac{\dot{a}}{a}{t^{\prime}}{\varphi^{\prime}}+\frac{2}{r}{r^{\prime}}{\varphi^{\prime}}+2\frac{\cos\theta}{\sin\theta}{\varphi^{\prime}}{\theta^{\prime}}$	$\displaystyle=0\,.$		(2.13)

Due to the six isometries of the FLRW spacetime, we can define conserved quantities associated with those.

There exist six Killing vectors for the metric (2.1) that satisfy the Killing equations $\nabla_{\mu}\xi_{\nu}+\nabla_{\nu}\xi_{\mu}=0$. Explicitly, these are given by

	$\displaystyle\xi_{(1),\mu}$	$\displaystyle=$	$\displaystyle a^{2}\left(0,\frac{1}{\sqrt{1-kr^{2}}}\sin{\theta}\cos{\varphi},r\sqrt{1-kr^{2}}\cos{\theta}\cos{\varphi},-r\sqrt{1-kr^{2}}\sin\theta\sin{\varphi}\right),$
	$\displaystyle\xi_{(2),\mu}$	$\displaystyle=$	$\displaystyle a^{2}\left(0,\frac{1}{\sqrt{1-kr^{2}}}\sin{\theta}\sin{\varphi},r\sqrt{1-kr^{2}}\cos{\theta}\sin{\varphi},r\sqrt{1-kr^{2}}\sin\theta\cos{\varphi}\right),$
	$\displaystyle\xi_{(3),\mu}$	$\displaystyle=$	$\displaystyle a^{2}\left(0,\frac{1}{\sqrt{1-kr^{2}}}\cos{\theta},-r\sqrt{1-kr^{2}}\sin{\theta},0\right),$
	$\displaystyle\xi_{(4),\mu}$	$\displaystyle=$	$\displaystyle a^{2}\left(0,0,r^{2}\sin\varphi,r^{2}\cos\theta\sin\theta\cos\varphi\right),$
	$\displaystyle\xi_{(5),\mu}$	$\displaystyle=$	$\displaystyle a^{2}\left(0,0,-r^{2}\cos\varphi,r^{2}\cos\theta\sin\theta\sin\varphi\right),$
	$\displaystyle\xi_{(6),\mu}$	$\displaystyle=$	$\displaystyle a^{2}\left(0,0,0,-r^{2}\sin^{2}\theta\right).$

In Euclidean space, these correspond to three translations and three rotations and are thus related to the conserved linear and angular momentum, respectively. The existence of Killing vectors implies the existence of conserved quantities $\mathcal{C}_{(i)}:=\xi_{(i),\mu}x^{\prime\mu}(\lambda)$ which are constant along the geodesics, that is,

	$\displaystyle\frac{a^{2}}{\sqrt{1-kr^{2}}}\sin{\theta}\cos{\varphi}r^{\prime}+a^{2}r\sqrt{1-kr^{2}}\cos{\theta}\cos{\varphi}\theta^{\prime}-a^{2}r\sqrt{1-kr^{2}}\sin\theta\sin{\varphi}\varphi^{\prime}$	$\displaystyle=p_{x}\,,$		(2.14)
	$\displaystyle\frac{a^{2}}{\sqrt{1-kr^{2}}}\sin{\theta}\sin{\varphi}r^{\prime}+a^{2}r\sqrt{1-kr^{2}}\cos{\theta}\sin{\varphi}\theta^{\prime}+a^{2}r\sqrt{1-kr^{2}}\sin\theta\cos{\varphi}\varphi^{\prime}$	$\displaystyle=p_{y}\,,$		(2.15)
	$\displaystyle\frac{a^{2}}{\sqrt{1-kr^{2}}}\cos\theta r^{\prime}-a^{2}r\sqrt{1-kr^{2}}\sin\theta\theta^{\prime}$	$\displaystyle=p_{z}\,,$		(2.16)
	$\displaystyle a^{2}r^{2}\sin\varphi\theta^{\prime}+a^{2}r^{2}\sin^{2}\theta\cot\theta\cos\varphi\varphi^{\prime}$	$\displaystyle=L_{x}\,,$		(2.17)
	$\displaystyle-a^{2}r^{2}\cos{\varphi}\theta^{\prime}+a^{2}r^{2}\sin^{2}\theta\cot\theta\sin\varphi\varphi^{\prime}$	$\displaystyle=L_{y}\,,$		(2.18)
	$\displaystyle-a^{2}r^{2}\sin^{2}{\theta}\varphi^{\prime}$	$\displaystyle=L_{z}\,.$		(2.19)

However, these are not all independent, only $p_{z}$, $L_{z}$, $p^{2}=p_{x}^{2}+p_{y}^{2}+p_{z}^{2}$ and $L^{2}=L_{x}^{2}+L_{y}^{2}+L_{z}^{2}$ are independent. A direct calculation shows that $L^{2}$ and $p^{2}$ take the neat forms

$$L^{2}=a^{4}r^{4}\left(\theta^{\prime}\right)^{2}+\csc^{2}\theta L_{z}^{2}\,,\qquad p^{2}=\frac{a^{4}\left(r^{2}\left(1-kr^{2}\right)^{2}\left(\left(\theta^{\prime}\right)^{2}+\sin^{2}(\theta)\left(\varphi^{\prime}\right)^{2}\right)+\left(r^{\prime}\right)^{2}\right)}{1-kr^{2}}\,.$$

(2.20)

It is convenient to rewrite these equations as follows

	$\displaystyle(\theta^{\prime})^{2}$	$\displaystyle=\frac{L^{2}}{a^{4}r^{4}}-\frac{L_{z}^{2}}{a^{4}r^{4}\sin^{2}\theta}\,,$		(2.21)
	$\displaystyle(r^{\prime})^{2}$	$\displaystyle=\frac{1-kr^{2}}{a^{4}}\left(p^{2}+kL^{2}-\frac{L^{2}}{r^{2}}\right).$		(2.22)

In this form, it is clear that one is dealing with two coupled first order differential equations. Note that, from the $\theta$-equation (2.21), we obtain a restriction for the $\theta$-motion, that is,

$$\arcsin\left(\frac{L_{z}}{L}\right)\leq\theta\leq\pi-\arcsin\left(\frac{L_{z}}{L}\right).$$

(2.23)

When $L_{z}=L$, the motion ensues in the equatorial plane, $\theta=\pi/2$. Finally, the $t$-component of the geodesic equation becomes

$$t^{\prime\prime}+(kL^{2}+p^{2})\frac{\dot{a}}{a^{3}}=0\implies(t^{\prime})^{2}-\frac{(kL^{2}+p^{2})}{a^{2}}=\mathcal{L},$$

(2.24)

where $\mathcal{L}=0$ for massless test particles and $\mathcal{L}=-1$ for massive test particles, respectively. We will present solutions to the equations (2.19), (2.21), (2.22) and (2.24) in Section 3.

As we are primarily interested in the study of azimuthal geodesics in a closed universe ($k=1$), let us investigate the full geodesic equation assuming $r=r_{0}=\text{constant}$. This means $r^{\prime\prime}=r^{\prime}=0$ which together with (2.11) implies $(1-r_{0})=0$ or $r_{0}=1$ since the other possible solution $r_{0}=0$ or both $\theta$ and $\varphi$ constant are of no interest to us. Hence, $r_{0}=1$ is the only possible solution for such geodesics.

When $r_{0}=1$ is used in (2.20), it implies $p=0$ and, together with (2.14)–(2.16), yields $p_{x}=p_{y}=p_{z}=0$. Consequently, Eq. (2.22) is also identically satisfied.

Next, we assume $\theta=\theta_{0}=\text{constant}$. When this is put into (2.12), one immediately arrives at $\cos\theta_{0}\sin\theta_{0}{\varphi^{\prime}}^{2}=0$. Since solutions with constant $\varphi$ are of no interest to us, we are left with $\theta_{0}=\pi/2$. We neglect the solution $\theta_{0}=0$ as it corresponds to a pole. This yields $L_{x}=L_{y}=0$ by (2.17) and (2.18), and $L=L_{z}=-a^{2}\varphi^{\prime}$. This final relation is the starting point for the study of azimuthal geodesics.

Beginning with (2.21) and (2.19), one can eliminate the radial coordinate and the temporal coordinate to arrive at

$$\derivative{\theta}{\varphi}=\sin{\theta}\sqrt{\frac{L^{2}}{L_{z}^{2}}\sin^{2}\theta-1}\,,$$

(3.1)

which can be integrated to give

$$\cot^{2}\negmedspace\theta=\left(\frac{L^{2}}{L_{z}^{2}}-1\right)\sin^{2}(\varphi-\varphi_{0})\,,$$

(3.2)

where $\varphi_{0}$ is a constant of integration, fixed by the initial conditions. Equation (3.2) describes arbitrary great circles, and we note that for $L=L_{z}$ one finds $\theta=\pi/2$. For all $L\neq L_{z}$, the geodesic motion is three-dimensional with the angular momentum $L$ precessing around the $z$-axis. In fact, due to the spatial isotropy one can always choose $L=L_{z}$ and hence the azimuthal geodesics are always great circles.

By employing (2.22) and (2.21), we find

$$\left(\derivative{r}{\theta}\right)^{2}=-\frac{r^{2}\left(kr^{2}-1\right)\left(L^{2}\left(kr^{2}-1\right)+p^{2}r^{2}\right)}{L^{2}-L_{z}^{2}\csc^{2}(\theta)}\,,$$

(3.3)

which can be re-written as

$$\left(\derivative{r}{\theta}\right)^{2}=-r^{2}\left[\left(k^{2}L^{2}+kp^{2}\right)r^{4}+\left(-2kL^{2}-p^{2}\right)r^{2}+L^{2}\right]\frac{1}{L^{2}}\left(1-\frac{1}{\left(L/L_{z}\right)^{2}\sin^{2}(\theta)}\right)^{-1}.$$

(3.4)

Now, (3.4) can be solved by separation of variables, that is,

$$\operatorname{arctanh}\left(\sqrt{\frac{p^{2}}{L^{2}}\frac{r^{2}}{\left(kr^{2}-1\right)}+1}\right)=\left[\arcsin\left(\frac{\cos{\theta}}{\sqrt{1-\left(L_{z}/L\right)^{2}}}\right)-\arcsin\left(\frac{\cos{\theta_{0}}}{\sqrt{1-\left(L_{z}/L\right)^{2}}}\right)\right],$$

(3.5)

where $\theta_{0}$ is fixed by the initial conditions and depends, in fact, on the orbit’s radius. If the orbit has a finite maximal radius (which is, for example, always true for $k=1$), it is convenient to choose $\theta_{0}=\pi-\arcsin(L/L_{z})$ so that (3.5) reads

$$\operatorname{arctanh}\left(\sqrt{\frac{p^{2}}{L^{2}}\frac{r^{2}}{\left(kr^{2}-1\right)}+1}\right)=\arcsin\left(\frac{\cos{\theta}}{\sqrt{1-\left(L_{z}/L\right)^{2}}}\right)+\frac{\pi}{2}\,.$$

(3.6)

If the orbit’s radius extends to infinity (which is possible, for example, when $k=0$ or $k=-1$), we choose $\theta_{0}=\pi/2$ such that (3.5) becomes

$$\operatorname{arctanh}\left(\sqrt{\frac{p^{2}}{L^{2}}\frac{r^{2}}{\left(kr^{2}-1\right)}+1}\right)=\arcsin\left(\frac{\cos{\theta}}{\sqrt{1-\left(L_{z}/L\right)^{2}}}\right).$$

(3.7)

We can now consider (2.22) and (2.24) to find a relation for $r(t)$ or, equivalently, $t(r)$. In particular, we have

$$\left(\derivative{r}{t}\right)^{2}=-\frac{\left(kr^{2}-1\right)\left(L^{2}\left(kr^{2}-1\right)+p^{2}r^{2}\right)}{r^{2}a^{2}\left(a^{2}\mathcal{L}+kL^{2}+p^{2}\right)}\,.$$

(3.8)

This equation is separable and we can find

$$\int\frac{r\,\mathrm{d}r}{\sqrt{-\left(kr^{2}-1\right)\left(L^{2}\left(kr^{2}-1\right)+p^{2}r^{2}\right)}}\\ =\frac{1}{\sqrt{k(kL^{2}+p^{2})}}\arctan(\sqrt{\frac{p^{2}}{\left(1-kr^{2}\right)\left(kL^{2}+p^{2}\right)}-1})+c_{1}$$

(3.9)

To find the integral with respect to $t$, we need to specify $a(t)$. We consider two general forms for $a(t)$ which are usually discussed in this context. First, consider a power law which arises e.g. for a matter or radiation dominated universe,

$$a(t)=\alpha t^{\mu}\,,$$

(3.10)

where $\alpha$ is an integration constant and $\mu\in\mathbb{R}^{+}$ is a dimensionless parameter. We therefore have

$$\int\frac{\mathrm{d}t}{a\sqrt{\left(a^{2}\mathcal{L}+kL^{2}+p^{2}\right)}}=-\frac{t^{1-\mu}\sqrt{kL^{2}+p^{2}+\alpha^{2}\mathcal{L}t^{2\mu}}\,_{2}F_{1}\left(1,\frac{1}{2\mu};\frac{\mu+1}{2\mu};-\frac{t^{2\mu}\alpha^{2}\mathcal{L}}{kL^{2}+p^{2}}\right)}{\alpha(\mu-1)\left(kL^{2}+p^{2}\right)}+c_{2}\,,$$

(3.11)

and we note that for massless test particles ($\mathcal{L}=0$), we find the relation

$$\frac{1}{\sqrt{k}}\arctan(\sqrt{\frac{p^{2}}{\left(1-kr^{2}\right)\left(kL^{2}+p^{2}\right)}-1})=\frac{t^{1-\mu}}{(\alpha-\alpha\mu)}+c_{3}.$$

(3.12)

Here $c_{2}$ and $c_{3}$ denote constants of integration.

Lastly, we consider an exponential form, e.g. for a universe to be dominated by a cosmological constant $\lambda>0$,

$$a(t)=\alpha\exp\left(\lambda t\right),$$

(3.13)

where $\alpha$ denotes an integration constant. Then we have

$$\int\frac{\mathrm{d}t}{a\sqrt{\left(a^{2}\mathcal{L}+kL^{2}+p^{2}\right)}}=-\frac{\mathrm{e}^{-\lambda t}\sqrt{kL^{2}+p^{2}+\alpha^{2}\mathcal{L}\mathrm{e}^{2\lambda t}}}{\alpha k\lambda L^{2}+\alpha\lambda p^{2}}+c_{2}\,.$$

(3.14)

We once again note that for $\mathcal{L}=0$, i.e.  for a massless test particle, we obtain the relation

$$\frac{1}{\sqrt{k}}\arctan(\sqrt{\frac{p^{2}}{\left(1-kr^{2}\right)\left(kL^{2}+p^{2}\right)}-1})=-\frac{\mathrm{e}^{-\lambda t}}{\alpha\lambda}+c_{3}\,.$$

(3.15)

In principle, this can be rewritten to find an explicit expression for the radial coordinate $r$ in terms of cosmological time $t$. Next, we will present a comprehensive study of the azimuthal geodesics.

As we can see from (2.5), the energy-momentum density evolves differently according to the equation of state parameter. In the case of a universe filled with two barotropic fluids, this means that one of the fluids may prevail for some time in the universe’s evolution. This approach allows for a more comprehensive understanding of the dynamics of the universe. For a review on multifluid cosmologies see e.g. [36]. Azimuthal geodesics for light-like particles in single-fluid cosmologies have been discussed previously [51]. Here, we return to the case of massless particles and now extend the approach to the two-fluid case.

For convenience, we use the rescaling

$$a\mapsto\left(\beta_{1}\right)^{1/\gamma_{1}}a\,,\quad\text{and}\quad t\mapsto\left(\beta_{1}\right)^{1/\gamma_{1}}t\,,$$

(4.4)

such that Eq. (2.6) can be re-written in terms of $\xi:=\beta_{2}/(\beta_{1})^{\gamma_{2}/\gamma_{1}}$, and reads

$$\left(\dot{a}\right)^{2}=a^{-\gamma_{1}}+\xi a^{-\gamma_{2}}-1\,.$$

(4.5)

Moreover, since $\mathcal{L}=0$ is for null geodesics, Eq. (4.3) becomes

$\displaystyle\Delta\varphi=2\int\limits_{0}^{a_{\rm max}}\left[a^{2}\left(a^{-\gamma_{1}}+\xi a^{-\gamma_{2}}-1\right)\right]^{-1/2}\mathrm{d}a\,.$

(4.6)

This integral can be given in terms of elementary functions only for specific choices of $\gamma_{1}$ and $\gamma_{2}$. This is related to the Theorem of Chebyshev, see Appendix A.1.

We will now discuss the azimuthal geodesics of a massless particle in a universe filled by matter, radiation, or stiff matter, together with a cosmological constant. To this end, we note that, for a cosmological constant $\Lambda$, we can use (2.6) where we set $\beta_{2}\equiv\Lambda/3$ and $\gamma_{2}=-2$, or, equivalently, $w_{2}=-1$, effectively treating the cosmological constant as a fluid. In the following, we also set $\beta_{1}=\beta$, $\gamma_{1}=\gamma$, and, thus, $\xi=\Lambda\beta^{2/\gamma}/3$. Note that $\xi\gtrless 0$ for $\Lambda\gtrless 0$.

Before proceeding with our analysis, we seek the condition that the solution $a_{\rm max}$ has to satisfy in order for it to be a positive real maximum. Equations (2.6) and (2.3) read

	$\displaystyle\dot{a}^{2}$	$\displaystyle=a^{-\gamma}+\xi a^{2}-1\,,$		(4.18)
	$\displaystyle a\ddot{a}$	$\displaystyle=-\frac{1}{2}\gamma a^{-\gamma}+\xi a^{2}=-\frac{1}{2}\gamma a^{-\gamma}+\xi a^{2}\,.$		(4.19)