Exact Geodesic Distances in FLRW Spacetimes

TL;DR

This paper develops methods to compute exact geodesic distances in flat FLRW spacetimes, providing closed-form solutions for simple cases and a fast numerical approach for complex models relevant to cosmology.

Contribution

It introduces a general framework for calculating geodesic distances in flat FLRW spacetimes, including closed-form solutions and a numerical method for complex cases.

Findings

Closed-form solutions for geodesic distances in flat spacetimes with dark energy or matter.

A fast numerical method for geodesics in mixed dark energy and matter spacetimes.

A general approach to determine geodesic connectedness based on the scale factor.

Abstract

Geodesics are used in a wide array of applications in cosmology and astrophysics. However, it is not a trivial task to efficiently calculate exact geodesic distances in an arbitrary spacetime. We show that in spatially flat (3+1)-dimensional Friedmann-Lemaitre-Robertson-Walker (FLRW) spacetimes, it is possible to integrate the second-order geodesic differential equations, and derive a general method for finding both timelike and spacelike distances given initial-value or boundary-value constraints. In flat spacetimes with either dark energy or matter, whether dust, radiation, or a stiff fluid, we find an exact closed-form solution for geodesic distances. In spacetimes with a mixture of dark energy and matter, including spacetimes used to model our physical universe, there exists no closed-form solution, but we provide a fast numerical method to compute geodesics. A general method is also described for determining the geodesic connectedness of an FLRW manifold, provided only its scale factor.

Equations

$$d{s}^2=-d{t}^2+a(t{)}^{2}d{\Sigma }^2,$$

$${\left(\frac{\dot{a}}{a}\right)}^2=\frac{\Lambda }{3}+\frac{c}{a^{3g}},$$

$$a_{\Lambda }(t)$$

$$a_D(t)$$

$$a_R(t)$$

$$a_S(t)$$

$$a_{\Lambda D}(t)$$

$$a_{\Lambda R}(t)$$

$$a_{\Lambda S}(t)$$

$$-z_0^2+z_1^2+z_2^2+z_3^2+z_4^2={\lambda }^2,$$

$$\begin{array}{rl}\frac{{\lambda }^2+s^2}{2\eta }& \mapsto z_0,\\ \frac{{\lambda }^2-s^2}{2\eta }& \mapsto z_1,\\ \frac{\lambda }{\eta }r\cos \theta & \mapsto z_2,\\ \frac{\lambda }{\eta }r\sin \theta \cos \varphi & \mapsto z_3,\\ \frac{\lambda }{\eta }r\sin \theta \sin \varphi & \mapsto z_4,\end{array}$$

$$\eta (t)={\int }^t\frac{d{t}^{\prime }}{a(t^{\prime })}.$$

$$d(x,y)=\{\begin{array}{ll}\lambda \mathrm{arccosh}\frac{⟨x,y⟩}{{\lambda }^2}& \text{if }x-y\text{ is timelike,}\\ 0& \text{if }x-y\text{ is lightlike,}\\ \infty & \text{if }⟨x,y⟩\le -{\lambda }^2\text{ and }x\ne -y,\\ \lambda \arccos \frac{⟨x,y⟩}{{\lambda }^2}& \text{otherwise.}\end{array}$$

$$\delta \int ds=0,$$

$$\delta {\int }_{{\sigma }_1}^{{\sigma }_2}\sqrt{g_{\mu \nu }\frac{\partial x^{\mu }}{\partial \sigma }\frac{\partial x^{\nu }}{\partial \sigma }}d\sigma =0.$$

$${\nabla }_X\frac{\partial x^{\mu }}{\partial \sigma }=\frac{{\partial }^2{x}^{\mu }}{\partial {\sigma }^2}+{\Gamma }_{\rho \tau }^{\mu }\frac{\partial x^{\rho }}{\partial \sigma }\frac{\partial x^{\tau }}{\partial \sigma }=\gamma \left(\sigma \right)\frac{\partial x^{\mu }}{\partial \sigma },$$

$${\Gamma }_{\rho \tau }^{\mu }=\frac{1}{2}{g}^{\mu \nu }\left(\frac{\partial g_{\nu \rho }}{\partial x^{\tau }}+\frac{\partial g_{\nu \tau }}{\partial x^{\rho }}-\frac{\partial g_{\rho \tau }}{\partial x^{\nu }}\right),$$

$$\frac{{\partial }^{2}t}{\partial {\sigma }^2}+a\frac{da}{dt}{h}_{ij}\frac{\partial x^i}{\partial \sigma }\frac{\partial x^j}{\partial \sigma }$$

$$\frac{{\partial }^2{x}^i}{\partial {\sigma }^2}+\frac{2}{a}\frac{da}{dt}\frac{\partial t}{\partial \sigma }\frac{\partial x^i}{\partial \sigma }+{\Gamma }_{jk}^i\frac{\partial x^j}{\partial \sigma }\frac{\partial x^k}{\partial \sigma }$$

$$\omega$$

$${\left(\frac{\partial \omega }{\partial \sigma }\right)}^2$$

$$\begin{array}{rl}\frac{\partial \omega }{\partial \sigma }\frac{{\partial }^2\omega }{\partial {\sigma }^2}& -\frac{1}{2}\frac{\partial h_{ij}}{\partial x^k}\frac{\partial x^i}{\partial \sigma }\frac{\partial x^j}{\partial \sigma }\frac{\partial x^k}{\partial \sigma }+\frac{2}{a}\frac{da}{dt}\frac{\partial t}{\partial \sigma }{\left(\frac{\partial \omega }{\partial \sigma }\right)}^2\\ & +h_{ij}{\Gamma }_{kl}^i\frac{\partial x^j}{\partial \sigma }\frac{\partial x^k}{\partial \sigma }\frac{\partial x^l}{\partial \sigma }=\gamma {\left(\frac{\partial \omega }{\partial \sigma }\right)}^2.\end{array}$$

$$\frac{{\partial }^{2}t}{\partial {\sigma }^2}+a\frac{da}{dt}{\left(\frac{\partial \omega }{\partial \sigma }\right)}^2$$

$$\frac{{\partial }^2\omega }{\partial {\sigma }^2}+\frac{2}{a}\frac{da}{dt}\frac{\partial t}{\partial \sigma }\frac{\partial \omega }{\partial \sigma }$$

$$\frac{{\partial }^{2}t}{\partial {\omega }^2}+\frac{da}{dt}\left(a-\frac{2}{a}{\left(\frac{\partial t}{\partial \omega }\right)}^2\right)=0.$$

$$\frac{{\partial }^{2}t}{\partial {\omega }^2}=-a\frac{da}{dt}.$$

$$0$$

$$\frac{\partial \omega }{\partial t}$$

$$D\left(\sigma \right)$$

$$d$$

$$D\left(\omega \right)$$

$$D\left(t\right)$$

$$D\left(t\right)=\sqrt{\frac{-\mu a^2\left(t\right)}{1+\mu a^2\left(t\right)}}.$$

$$\mu =a_0^{-2}\left(v_0^{-2}{a}_0^{-2}-1\right),$$

$$a(t_c)=\pm \sqrt{-{\mu }^{-1}}.$$

$${\omega }_c={\int }_{t_1}^{t_2}\frac{1}{a(t)}{\left(1-{\left(\frac{a(t)}{a(t_2)}\right)}^2\right)}^{-1/2}dt.$$

Full text

I Introduction

Cosmic microwave background experiments such as COBE G.F. Smoot et. al. (1992), WMAP G. Hinshaw et. al. (2013), and Planck P.A.R. Ade et. al. (2016) (Planck Collaboration) provide evidence for both early time cosmic inflation Guth (1981); Linde (1982) and late time acceleration S. Perlmutter et. al. (1998); A.G. Riess et. al. (1998), with interesting dynamics in between explaining many features of the universe, many of which are remarkably accurately predicted by the $\Lambda$CDM model Peebles (1982); Turner et al. (1984); Blumenthal et al. (1984); Davis et al. (1985). These and other experiments in recent decades have demonstrated that, to a high degree of precision, at large scales the visible universe is spatially homogeneous, isotropic, and flat, i.e., that its spacetime is described by the Friedmann-Lemaître-Robertson-Walker (FLRW) metric. FLRW spacetimes are therefore of particular interest in modern cosmology.

Here we develop a method for the exact calculation of the geodesic distance between any given pair of events in any flat FLRW spacetime. Geodesics and geodesic distances naturally arise in a wide variety of investigations not only in cosmology, but also in astrophysics and quantum gravity, with topics ranging from the horizon and dark energy problems, to gravitational lensing, to evaluating the observational signatures of cosmic bubble collisions and modified gravity theories, to the AdS/CFT correspondence Ellis and van Elst (1999); Grøn and Elgarøy (2007); Albareti et al. (2012); Demianski et al. (2003); Hirata and Seljak (2005); Bikwa et al. (2012); Doplicher et al. (2013); Melia (2013a); Bhattacharya and Tomaras (2017); Pyne and Birkinshaw (1996); Park (2008); Sereno (2009); Mukohyama (2009a, b); Traschen and Eardley (1986); Futamase and Sasaki (1989); Cooperstock et al. (1998); Caldwell and Langlois (2001); Dappiaggi et al. (2008); Kaloper et al. (2010); Melia (2013b); Bahrami (2017); Koyama and Soda (2001); Dong et al. (2012a, b); Minton and Sahakian (2008); Wainwright et al. (2014); Hagala et al. (2016). Closed-form solutions of the geodesic equations are also quite useful in validating a particular FLRW model by investigating how curvature, quintessence, local shear terms, etc., affect observational data. These solutions are of perhaps the greatest and most direct utility in large-scale N-body simulations, e.g., studying the large scale structure formation, which can benefit greatly from using such solutions by avoiding the costly numerical integration of the geodesic differential equations Adamek et al. (2016); Koksbang and Hannestad (2015); Bibiano and Croton (2017).

For a general spacetime, solving the geodesic equations exactly for given initial-value or boundary-value constraints is intractable, although it may be possible in some cases. For example, in $(3+1)$-dimensional de Sitter space, which represents a spacetime with only dark energy and is a maximally symmetric solution to Einstein’s equations, it turns out to be rather simple to study geodesics by embedding the manifold into flat $(4+1)$-dimensional Minkowski space $\mathbb{M}^{5}$. This construction was originally realized by de Sitter himself de Sitter (1917), and was later studied by Schrödinger Schrödinger (1956). In Sec. II.1 we review how geodesics may be found using the unique geometric properties of this manifold.

However, it is not so easy to explicitly calculate geodesic distances in other FLRW spacetimes except under certain assumptions. One approach would be to follow de Sitter’s philosophy by finding an embedding into a higher-dimensional manifold. Such an embedding always exists due to the Campbell-Magaard theorem, which states that any analytic $n$-dimensional Riemannian manifold may be locally embedded into an $(n+1)$-dimensional Ricci-flat space Campbell (1926); Magaard (1963), combined with a theorem due to A. Friedman extending the result to pseudo-Riemannian manifolds Friedman (1961, 1965). In fact, the embedding map is given explicitly by J. Rosen in Rosen (1965). However, it has since been shown that the metric in the embedding space is block diagonal with respect to the embedded surface, i.e., when the geodesic is constrained to the $(3+1)$-dimensional subspace we regain the original $(3+1)$-dimensional geodesic differential equations and we learn nothing new Romero et al. (1996).

Instead, in Sec. III we solve directly the geodesic differential equations for a general FLRW spacetime in terms of the scale factor and a set of initial-value or boundary-value constraints. The final geodesic distance can be written as an integral which is a function of the boundary conditions and one extra constant $\mu$, defined by a transcendental integral equation. This constant proves to be useful in a number of ways: it tells us if a manifold is geodesically connected provided only the scale factor. The solution of the integral equation defining this constant exists only if a geodesic exists for a given set of boundary conditions, and it helps one to find the geodesic distance, if such a geodesic exists.

We then give some examples in Sec. IV to show that for many scale factors of interest, we can find a closed-form solution. In cases where no closed-form solution exists, we can still transform the problem into one which is suitable for fast numerical integration. Concluding remarks are in Sec. V.

II Review of FLRW Spacetimes and de Sitter Embeddings

Friedmann-Lemaître-Robertson-Walker (FLRW) spacetimes are spatially homogeneous and isotropic $(3+1)$-dimensional Lorentzian manifolds which are solutions to Einstein’s equations Griffiths and Podolský (2009). These manifolds have a metric $g_{\mu\nu}$ with $\mu$, $\nu\in\{0,1,2,3\}$ that, when diagonalized in a given coordinate system, gives an invariant interval $ds^{2}=g_{\mu\nu}\,dx^{\mu}\,dx^{\nu}$ of the form

[TABLE]

where $a(t)$ is the scale factor, which describes how space expands with time $t$, and $d\Sigma$ is the spatial metric given by $d\Sigma^{2}=dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta\,d\phi^{2})$ for flat space in spherical coordinates $(r,\theta,\phi)$ that we use hereafter. The scale factor is found by solving Friedmann’s equation, the differential equation given by the $\mu=\nu=0$ component of Einstein’s equations:

[TABLE]

where $\Lambda$ is the cosmological constant, $g$ parametrizes the type of matter within the spacetime, $c$ is a constant proportional to the matter density, and we have assumed spatial flatness in our choice of $d\Sigma$. The scale factors for manifolds which represent spacetimes with dark energy ($\Lambda$), dust ($D$), radiation ($R$), a stiff fluid ($S$),111Stiff fluids are exotic forms of matter which have a speed of sound equal to the speed of light. They have been studied in a variety of models of the early universe, including kination fields, self-interacting (warm) dark matter, and Hor̂ava-Lifshitz cosmologies Dutta and Scherrer (2010). or some combination (e.g., $\Lambda D$ for dark energy and dust matter) are given by Griffiths and Podolský (2009)

[TABLE]

where $\lambda$ and $\alpha\equiv(c\lambda^{2})^{1/3}$ are the temporal and spatial scale-setting parameters. In manifolds with dark energy, i.e., $\Lambda>0$, $\lambda\equiv\sqrt{3/\Lambda}$.

II.1 de Sitter Spacetime

The de Sitter spacetime is one of the first and best studied spacetimes: de Sitter himself recognized that the $(3+1)$-dimensional manifold d$\mathbb{S}^{4}$ can be visualized as a single-sheet hyperboloid embedded in $\mathbb{M}^{5}$, defined by

[TABLE]

where $\lambda$ is the pseudo-radius of the hyperboloid de Sitter (1917). The injection $\chi\,:\,\mathrm{d}\mathbb{S}^{4}\hookrightarrow\mathbb{M}^{5}\,,\,\chi(x)\mapsto z$ is

[TABLE]

where $s^{2}\equiv r^{2}-\eta^{2}$, and the conformal time $\eta$ is defined as

[TABLE]

This embedding is a particular instance of the fact that any analytic $n$-dimensional pseudo-Riemannian manifold may be isometrically embedded into (at most) a $(n(n+1)/2)$-dimensional pseudo-Euclidean manifold (i.e., a flat metric with arbitrary non-Riemannian signature) Friedman (1961, 1965). The minimal $(n+1)$-dimensional embedding is most easily obtained using group theory by recognizing that the Lorentz group SO(1,3) is a stable subgroup of the de Sitter group dS(1,4) while the pseudo-orthogonal group SO(1,4) acts as its group of motions, i.e., dS(1,4) = SO(1,4)/SO(1,3), thereby indicating the minimal embedding is into the $\mathbb{M}^{5}$ space Aldrovandi and Pereira (1995).

If a spacelike geodesic extends far enough, there will exist an extremum, identified as $P_{3}$ in (c). As a result, it is simplest to use the spatial distance $\omega$ to parametrize these geodesics, though time can be used as well so long as those geodesics with turning points are broken into two parts at the point $P_{3}$.

Furthermore, it can be shown that geodesics on a de Sitter manifold follow the lines defined by the intersection of the hyperboloid with a hyperplane in $\mathbb{M}^{5}$ containing the origin and both endpoints of the geodesic Asmus (2009). An illustration of both timelike and spacelike geodesics constructed this way in d$\mathbb{S}^{2}$ embedded in $\mathbb{M}^{3}$ can be found in Fig. 1. This construction implies that the geodesic distance $d(x,y)$ in d$\mathbb{S}^{4}$ between two points $x$ and $y$ can be found using their inner product $\langle x,y\rangle=-x_{0}y_{0}+x_{1}y_{1}+x_{2}y_{2}+x_{3}y_{3}+x_{4}y_{4}$ in $\mathbb{M}^{5}$ via the following expression:

[TABLE]

While there are many ways to find geodesic distances on a de Sitter manifold, this is perhaps the simplest one.

III The Geodesic Equations in Four Dimensions

While de Sitter symmetries cannot be exploited in a general FLRW spacetime, it is still possible to solve the geodesic equations. A geodesic is defined in general by the variational equation

[TABLE]

which, if parametrized by parameter $\sigma$ ranging between two points $\sigma_{1}$ and $\sigma_{2}\,,$ becomes

[TABLE]

The corresponding Euler-Lagrange equations obtained via the variational principle yield the well-known geodesic differential equations:

[TABLE]

for the geodesic path $x^{\mu}(\sigma)$ with some as yet unknown function $\gamma(\sigma)$, where $\Gamma_{\rho\tau}^{\mu}$ are the Christoffel symbols defined by

[TABLE]

and $\nabla_{X}$ indicates the covariant derivative with respect to the tangent vector field $X$ O’Neill (1983). If the parameter $\sigma$ is affine, then $\gamma(\sigma)=0$. To solve a particular problem with constraints, we must use both (9) and (10).

III.1 The Differential Form of the Geodesic Equations

If only the non-zero Christoffel symbols are kept, then (10) can be broken into two differential equations written in terms of the scale factor:

[TABLE]

where $h_{ij}$ is the first fundamental form, i.e., the induced metric on a constant-time hypersurface, and the Latin indices are restricted to $\{1,2,3\}$.

To solve these, consider the spatial (Euclidean) distance $\omega$ between two points $\sigma_{1}$ and $\sigma_{2}$:

[TABLE]

This relation implies the spatial coordinates obey a geodesic equation with respect to the induced metric $h_{ij}$. Now, (12a) may be written in terms of $\omega$ using (13). The transformation needed for (12b) is found by multiplying by $h_{ij}(\partial x^{j}/\partial\sigma)$ and substituting the derivative of (13) with respect to $\omega$:

[TABLE]

The second and fourth terms cancel by symmetry and, supposing $(\partial\omega/\partial\sigma)\neq 0$, the pair of equations (12) may be written as

[TABLE]

We now proceed by parametrizing the geodesic by the Euclidean spatial distance, i.e., $\sigma\equiv\omega$. This yields $\partial\omega/\partial\sigma=1$, $\partial^{2}\omega/\partial\sigma^{2}=0$, and then (15b) gives $\gamma=(2/a)(da/dt)(\partial t/\partial\omega)$. Using these new relations, (15a) can be written as

[TABLE]

While neither the spatial distance $\omega$ nor time $t$ are affine parameters along all Lorentzian geodesics, the results will not be affected, since the differential equations no longer refer to $\gamma$. We can see that if $\partial t/\partial\omega=0$ then the second derivative of $t$ is always negative for $t>0$ and positive for $t<0$, since $da/dt>0$ for expanding spacetimes:

[TABLE]

If there exists a critical point exactly at $t=0$, it is a saddle point. From these facts, we conclude that any extremum found along a geodesic on a Friedmann-Lemaître-Robertson-Walker manifold is a local maximum in $t>0$ and a local minimum in $t<0$ with respect to $\omega$.222This statement is true under the assumption that the scale factor is a well-behaved monotonic function. If this condition does not hold, the following analysis must be reinspected. An example of such a curve with an extremum is shown in Fig. 1(c).

The second-order equation (16) may be simplified by multiplying by $2a^{-4}(\partial t/\partial\omega)$ and integrating by parts to get a non-linear first-order differential equation and a constant of integration $\mu$:

[TABLE]

the right hand side of which is hereafter referred to as the geodesic kernel $G(t;\mu)$. We may neglect the sign by noting that the spatial distance $\omega$ should always be an increasing function of $t$, so that any integration of the geodesic kernel should be always be performed from past to future times. It will prove necessary to know the value of $\mu$ to find the final value of the geodesic length between two events.

III.2 The Integral Form of the Geodesic Equations

To find the geodesic distance between a given pair of points/events, we need to use (19) in conjunction with the integral form of the geodesic equation, given in (9). We begin by defining the integrand in (9) as the distance kernel $D(\sigma)$:

[TABLE]

The invariant interval (1) tells us that $D^{2}$ is negative for timelike-separated pairs and positive for spacelike-separated ones, assuming $\sigma$ is monotonically increasing along the geodesic. Therefore, we always take the absolute value of $D^{2}$ so that the distance kernel is real-valued, while keeping in mind which type of geodesic we are discussing.

Depending on the particular scale factor and boundary values, we might sometimes parametrize the system using the spatial distance and other times using time. If we parametrize the geodesic with the spatial distance we find

[TABLE]

where the function $t(\omega)$ in the former equation is the inverted solution $\omega(t)$ to the differential equation (19). Since the distance kernel is a function of the geodesic kernel, we will need to know the value $\mu$ associated with a particular set of constraints.

If we insert (19) into (23), we can can see what values the constant $\mu$ can take:

[TABLE]

If $D^{2}<0$ for timelike intervals, then $\mu>0$. If $\mu=0$, we obtain a lightlike geodesic, since the distance kernel becomes zero. Hence, spacelike intervals correspond to $-a^{-2}(t)<\mu<0$. We do not consider $\mu<-a^{-2}(t)$ because this corresponds to an imaginary $\partial\omega/\partial t$, which we consider non-physical.

III.3 Geodesic Constraints and Critical Points

We would like to find geodesics for both initial-value and boundary-value problems. If we have Cauchy boundary conditions, i.e., the initial position and velocity vector are known, then finding $\mu$ is simple: since the left hand side of (19) is just the speed $v_{0}\equiv|v^{i}(t_{0})|$, where $v^{i}$ is the velocity vector defined by our initial conditions, we have

[TABLE]

where $a_{0}\equiv a(t_{0})$. This allows for simple solutions to cases with Cauchy boundary conditions.

However, if we have Dirichlet boundary conditions, i.e., the initial and final positions are known, then we must integrate (19) instead. The bounds of such an integral need to be carefully considered: if we have a spacelike geodesic which starts and ends at the same time, for instance, then it is not obvious how to integrate the geodesic kernel. In fact, we face an issue with the boundaries whenever we have geodesics with turning points. This feature occurs whenever $\partial t/\partial\omega=0$, i.e.,

[TABLE]

Since in this case $\mu<0$, we see that turning points only occur for spacelike geodesics. Furthermore, since all of the scale factors given by (3) are monotonic, this situation occurs in such spacetime only at a single point along a geodesic, if at all, identified as $P_{3}$ in Fig. 1(c). Specifically, if $t_{1},t_{2},t_{3}$ respectively correspond to the times at points $P_{1},P_{2},P_{3}$, then the integral of the geodesic kernel is found by integrating from $t_{3}$ to $t_{1}$ as well as from $t_{3}$ to $t_{2}$, since time is not monotonic along the geodesic. If no such turning point $P_{3}$ exists along the geodesic, a single integral from $t_{1}$ to $t_{2}$ may be performed. The integral of the distance kernel should be performed in the same way for the same reasons.

To determine if a turning point exists along a spacelike geodesic, we begin by noting that there is a corresponding critical spatial distance $\omega_{c}$ which corresponds to the critical time defined in (26). If we suppose $t_{2}>t_{1}>0$, then the geodesic kernel is maximized when $\mu=\mu_{c}\equiv-a^{-2}(t_{2})$, i.e., when $\mu$ attains its minimum value. This is the minimum value of $\mu$ along the geodesic, since $a(t)$ is monotonically increasing. The critical spatial distance is defined by this $\mu_{c}$ and is given by

[TABLE]

Since $\mu_{c}$ maximizes the geodesic kernel, it is impossible for a spacelike-separated pair to be spatially farther apart without their geodesic having a turning point. We then conclude that if $\omega<\omega_{c}$ for a particular pair of spacelike-separated points, then the geodesic is of the form shown in Fig. 1(b), and if $\omega>\omega_{c}$ it is of the form shown in Fig. 1(c). In other words, if the geodesic is of the latter type, then the solution to (19) is

[TABLE]

again supposing $t_{2}>t_{1}>0$. The bounds on the integral are chosen this way due to the change of sign in the geodesic kernel on opposite sides of the critical point. If $0>t_{2}>t_{1}$ then the bounds on the integrals are reversed so that $\omega,d>0$.

III.4 Geodesic Connectedness

Certain FLRW manifolds are not spacelike-geodesically-connected, meaning not all pairs of spacelike-separated points are connected by a geodesic. For a given pair of times $t_{1},t_{2}$ there exists a maximum spatial separation $\omega_{m}$ past which the two points cannot be connected by a geodesic. To determine this maximum spatial distance $\omega_{m}$ for a particular pair of points, we use (28), this time taking the limit $\mu\to 0^{-}$. This limit describes a spacelike geodesic which is asymptotically becoming lightlike. If the critical time $t_{c}$ remains finite in this limit, the manifold is geodesically connected and $\omega_{m}=\infty$, whereas if it becomes infinite then $\omega_{m}$ remains finite, shown in detail in Fig. 2. The equation (28) in the limit $\mu\to 0^{-}$ is

[TABLE]

Comparing (30) to (6) we notice that $\omega_{m}$ is simply a combination of conformal times using the boundary points $t_{1}$ and $t_{2}$: $\omega_{m}\propto\eta_{c}\equiv\eta(t_{c})$, and so if $\eta_{c}$ is finite, then $\omega_{m}$ will be finite as well. Therefore, we conclude that a Friedmann-Lemaître-Robertson-Walker manifold is geodesically complete if

[TABLE]

using the appropriate $a(t)$ and $t(\eta)$ for the given manifold.

As an example, consider the de Sitter manifold:

[TABLE]

so that the limit maximum conformal time in terms of $\mu$ is

[TABLE]

Therefore, in the flat foliation, there exist pairs of points on the de Sitter manifold which cannot be connected by a geodesic. On the other hand, if we consider the Einstein-de Sitter manifold, which represents a spacetime with dust matter, the scale factor is proportional to $\eta^{2}$:

[TABLE]

so that every pair of points may be connected by a geodesic.

IV Examples

Here we apply the results above to calculate geodesics in the FLRW manifolds defined by each of the scale factors in (3), using two type of constraints: the Dirichlet and Cauchy boundary conditions. The former conditions specify two events or points in a given spacetime that can be either timelike or spacelike separated, as in Fig. 1. The latter conditions specify just one point and a vector of initial velocity. If the initial speed is below the speed of light, then the resulting geodesic is timelike, and corresponds to a possible world line of a massive particle. If the initial speed is above the speed of light, i.e., the initial tangent vector is spacelike, then the resulting geodesic is spacelike, and corresponds to a geodesic of a hypothetical superluminal particle. Even though tachyons may not exist, spacelike geodesics are well defined mathematically. The last example that we consider illustrates how to apply these techniques to find numerical values for geodesic distances in our physical universe.

IV.1 Dark Energy

Suppose we wish to find the geodesic distance using the Dirichlet boundary conditions $\{t_{1},\allowbreak t_{2},\allowbreak\omega\}$. The geodesic kernel in a flat de Sitter spacetime is

[TABLE]

where $\mu$ has absorbed a factor of $\lambda^{2}$ and we use $\eta\in[-1,0)$ so that $t\geq 0$. We can easily transform the kernel into a polynomial equation by using the conformal time:

[TABLE]

If the minimal value of $\mu$ is inserted into this kernel, the turning point $\omega_{c}$ can be found exactly:

[TABLE]

The geodesic kernel may now be integrated both above and below the turning point:

[TABLE]

The variable $\mu$ is then found by inverting one of these equations. Finally, substitution of the scale factor and numerical value $\mu$ into (23) gives the geodesic distance for a pair of coordinates defined by $\{\eta_{1},\allowbreak\eta_{2},\allowbreak\omega\}$:

[TABLE]

for spacelike-separated pairs. In Appendix Appendix: Equivalence of de Sitter Solutions we show that this solution is equivalent to the solution found using the embedding in Sec. II.1.

IV.2 Dust

In this example, let us suppose we have Cauchy boundary conditions and we want an expression for the geodesic distance in terms of spatial distance traveled $\omega$. First, knowing the values $(t_{0},r_{0},\theta_{0},\phi_{0})$ and $|v_{0}|$, we can find the parameter $\mu$ via (25). Because the manifold has a singularity at $t=0$, we assert $t_{0}\neq 0$ to avoid a nonsensical value for $\mu$. We proceed by parametrizing the geodesic equation by the spatial distance, following (22), so that the distance kernel for this spacetime is

[TABLE]

We use the geodesic kernel to find $t(\omega)$ directly, by solving (19) for $\omega(t)$ and inverting the solution. In the spacetime with dust matter and no cosmological constant the geodesic kernel is

[TABLE]

which, using the transformations $x\equiv(3t/2\lambda)^{1/3}$ and $\mu\to\alpha^{2}\mu$, becomes

[TABLE]

The value of $\omega$ where the turning point occurs is then

[TABLE]

where $x_{0}\equiv x(t_{0})$ and ${}_{2}F_{1}(a,b;c;z)$ is the Gauss hypergeometric function.

The final expression $\omega(t)$ still depends on the existence of a critical point along the geodesic. To demonstrate how piecewise solutions are found, hereafter we suppose we are studying a superluminal inertial object moving fast and long enough to take a geodesic with a turning point. The spatial distance $\omega(x;x_{0})$, with $x>x_{0}$ and $\mu<0$, which we know because the geodesic is spacelike, is

[TABLE]

afterward, where $x_{c}=(-\mu)^{-1/4}$, and $K(m)$ and $F(\phi|m)$ respectively are the complete and incomplete elliptic integrals of the first kind with parameter $m$. These expressions $\omega(x;x_{0},x_{c})$ are slightly different for $\mu>0$. Despite the apparent complexity of the above expressions, they are in fact easy to invert via the Jacobi elliptic functions. The distance for a geodesic with a turning point is

[TABLE]

giving the final result

[TABLE]

where we have used the auxiliary variables

[TABLE]

and the three functions in the last definition are the Jacobi elliptic functions with parameter $m$.

IV.3 Radiation

Here we suppose we have Cauchy boundary conditions, but the particle will take a timelike geodesic, i.e., $\mu>0$. Using the transformations $x\equiv\sqrt{2t/\lambda}$ and $\mu\to\alpha^{3/2}\mu$, we can write the geodesic kernel as

[TABLE]

If this kernel is integrated over $x$ to find the spatial distance $\omega(x)$, the result can be inverted to give

[TABLE]

where $\beta_{1}\equiv\alpha^{3/4}\mu^{1/2}/\lambda$ and $\beta_{2}\equiv\operatorname{arcsinh}(\mu^{1/2}x_{0})$. Since the geodesic distance is more easily found when we parametrize with the spatial distance $\omega$, we can write the distance kernel as

[TABLE]

and the geodesic distance as

[TABLE]

Typically, timelike geodesics are parametrized by time: since there exists a closed-form solution for $\omega(x(t))$ this expression can be substituted here, though it would needlessly add extra calculations. Therefore, in practice it is computationally simpler to use a spatial parametrization.

IV.4 Stiff Fluid

Suppose we have a spacetime containing a homogeneous stiff fluid, and we wish to find a timelike geodesic using Dirichlet boundary conditions. Using the transformation $x\equiv(3t/\lambda)^{1/3}$, we can write the geodesic kernel as

[TABLE]

where $\mu$ has absorbed a factor of $\alpha$. This kernel can easily be integrated to find

[TABLE]

The constant $\mu$ may be found provided the initial conditions $\{x_{0},x_{1},\omega\}$:

[TABLE]

Finally, if the geodesic is parametrized by $x(t)$ we arrive at

[TABLE]

IV.5 Dark Energy and Dust

None of the spacetimes with a mixture of dark energy and some form of matter have closed-form solutions for geodesics, because the scale factors are various powers of the hyperbolic sine function, so it becomes cumbersome to work with the geodesic and distance kernels. However, by using the right transformations, it is still possible to make the problem well-suited for fast numerical integration. In this example, we use the mixed dust and dark energy spacetime, following the same procedure as before; for other spacetimes with mixed contents the same method applies. This time, the geodesic kernel is

[TABLE]

Once again, the kernel can be written as a polynomial expression, this time using the square root of the scale factor as the transformation:

[TABLE]

There is no known closed-form solution to the integral of $G_{\Lambda D}$. The distance kernel is best represented as a function of $t$ to simplify numerical evaluations:

[TABLE]

There is no known closed-form solution to this kernel’s integral either, but it can be quickly computed numerically, since the hyperbolic term needs to be evaluated only once for each value of $t$. In general, the numeric evaluations of such integrals can be quite fast if the kernels take a polynomial form, and a Gauss-Kronrod quadrature can be used for numeric evaluation of these integrals.

IV.6 Dark Energy, Dust, and Radiation

Typically in cosmology one studies one particular era, whether the early inflationary phase, the radiation-dominated phase, the matter-dominated phase after recombination, or ultimately today’s period of accelerated expansion. Perhaps the most important spacetime which we have not looked at yet is the FLRW spacetime which most closely models our own physical universe, in its entirety. In this section we will show how to most efficiently find geodesics in our (FLRW $\Lambda$DR) universe.

Because the scale factor $a(t)$ is a smooth, monotonic, differentiable, and bijective function of time, it, instead of time $t$ or spatial distance $\omega$, can parametrize geodesics, so long as we remember to break up expressions when there exists a turning point in long spacelike geodesics. In what follows we will restrict the analysis to timelike geodesics for simplicity. To find spacelike geodesics, refer to the steps performed in Sec. IV.2. Using the scale-factor parametrization, the geodesic and distance kernels are

[TABLE]

As we saw in Sec. IV.5, integrands such as these produce no closed-form solutions, but they are easily evaluated numerically due to their polynomial form.

We now provide a simple example of computing an exact geodesic distance between a pair of events in our physical universe using these results. Suppose we are to measure the timelike geodesic distance between an event in the early universe, where $t_{1}=10^{11}$s, and another event near today, $t_{2}=4.3\times 10^{17}$s. Let the spatial distance of this geodesic be $\omega=4.1\times 10^{13}$km, roughly the distance to Alpha Centauri. Taking relevant experimental values from recent measurements E. Calabrese et. al. (2017), we find the Hubble constant is $H_{0}=100h$ km/s/Mpc, where $h=0.705$, and the density parameters are $\Omega_{\Lambda}=0.723$, $\Omega_{D}=0.277$, and $\Omega_{R}=9.29\times 10^{-5}$. The leading constant $\lambda$ in the above equations can be expressed as $\lambda=H_{0}^{-1}\Omega_{\Lambda}^{-1/2}$, thereby completing the set of all the relevant physical parameters used in (66) and (67). We then integrate the geodesic kernel (66), inserting the speed of light $c$ where needed, to numerically solve for the integration constant $\mu$, which we find to be $\mu=2.53\times 10^{23}$. Inserting this value into the distance kernel (67) and evaluating numerically gives a final geodesic distance of $d=2.22\times 10^{23}$km.

V Conclusion

By integrating the geodesic differential equations (10) we have shown for spacetimes with dark energy, dust, radiation, or a stiff fluid, that it is possible to find a closed-form solution for the geodesic distance provided either initial-value or boundary-value constraints. Furthermore, by studying the form of the first-order differential equation (19) we found that extrema along spacelike geodesic curves will always point away from the origin. This insight provides a better understanding of how to integrate the geodesic and distance kernels (19, 22, 23) for different types of boundary conditions. Moreover, our other important result in Sec. III.4 demonstrates how, using (6), (26) and (31), we are able to tell, using only the scale factor, whether or not all points on a flat FLRW manifold can be connected by a geodesic. This observation is particularly useful in numeric experiments and investigations that can study only a finite portion of a spatially flat manifold. Finally, in Sec. IV we provided several examples of how these results might be applied to some of the most well-studied FLRW manifolds, including the manifold describing our universe. While not all spacetimes have closed-form solutions for geodesics, it is still possible to reframe the problem in a way which may be solved efficiently using numerical methods in existing software libraries.

Acknowledgements.

We thank Cody Long, Aron Wall, and Michel Buck for useful discussions and suggestions. This work was supported by NSF grants No. CNS-1442999, CCF-1212778, and PHY-1620526, ARO grant No. W911NF-16-1-0391, and DARPA grant No. N66001-15-1-4064. Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of DARPA.

Appendix: Equivalence

of de Sitter Solutions

Here we show that the equations (7) and (42) are equal under certain assumptions. Let us refer to the former as $d_{1}$ and the latter as $d_{2}$. The conformal time in the de Sitter spacetime is $\eta(t)=-e^{-t/\lambda}$, with $\eta\in[-1,0)$ so that the cosmological time $t$ remains positive. Since the geodesic distance depends on the spatial distance, but not the individual spatial coordinates, we can assume without loss of generality that the initial point is located at the origin, $r=\theta=\phi=0$, and the second point is located at some distance $\omega$ from the origin, $r=\omega,\,\theta=\phi=0$. Further, to simplify the proof, suppose the initial point is at time $t=0$ ($\eta=-1$) and the second point at some $t=t_{0}>0$ ($\eta=\eta_{0}\in(-1,0)$). We are allowed to make these assumptions due to the spatial symmetries associated with the dS(1,3) group and the existence of a global timelike Killing vector in the flat foliation of the de Sitter manifold Podolský (1993). In addition, suppose the geodesic is timelike so that $\omega\in[0,\eta_{0}+1)\subseteq[0,1)$. This same method may be applied to spacelike geodesics.

Using these values, the embedding coordinates in $\mathbb{M}^{5}$ are $x=((1-\lambda^{2})/2,\,-(1+\lambda^{2})/2,\,0,\,0,\,0)$ and $y=((\lambda^{2}+\omega^{2}-\eta_{0}^{2})/2\eta_{0},\,(\lambda^{2}-\omega^{2}+\eta_{0}^{2})/2\eta_{0},\,\lambda\omega/\eta_{0},\,0,\,0)$. This equation gives a geodesic distance

[TABLE]

On the other hand, we can use the solution provided by (42) using the value of $\mu$ in (41):

[TABLE]

in the geodesic distance expression

[TABLE]

If we apply $\cosh(d/\lambda)$ to each of these expressions, and use the identities $\cosh(x-y)=\cosh x\cosh y-\sinh x\sinh y$ and $\cosh\operatorname{arcsinh}x=\sqrt{x^{2}+1}$, we may equate them to get

[TABLE]

Using (69) and some algebra, the right hand side may be simplified to give the result on the left hand side, thereby proving they are equal.