Constant acceleration over cosmological distances

I've never seen the cosmological version of the relativistic rocket in any textbook, but I think it's fairly straightforward to derive it from the standard cosmological equations.

Let's start with the FLRW metric, $$ \text{d}s^2 = c^2\text{d}t^2 - a(t)\,\text{d}\ell^2, $$ where $a(t)$ is the scale factor and $\text{d}\ell$ the infinitesimal co-moving distance. The Friedmann equations for the standard ΛCDM-model have the solution $$ H(a) = \frac{\dot{a}}{a} = H_0\sqrt{\Omega_{R,0}\,a^{-4} + \Omega_{M,0}\,a^{-3} + \Omega_{K,0}\,a^{-2} + \Omega_{\Lambda,0}}, $$ which expresses the Hubble parameter $H(a)$ as a function of the Hubble constant and the relative present-day radiation, matter, and dark energy densities. From $$ \dot{a} = \frac{\text{d}a}{\text{d}t} $$ we get $$ \text{d}t = \frac{\text{d}a}{\dot{a}} = \frac{\text{d}a}{a\,H(a)}, $$ so that $$ t(a) = \int_0^a\frac{\text{d}a'}{a'\,H(a')}, $$ which we can numerically invert to obtain $a(t)$ (see also this post).

Now, a rocket with velocity $v(t)$ will travel in a time $\text{d}t$ a proper distance $$ a(t)\,\text{d}\ell = v(t)\,\text{d}t, $$ so that the total co-moving distance travelled in a cosmic time interval $[t_0,t_1]$ is given by $$ D_\text{c} = \int_{\ell_0}^{\ell_1}\text{d}\ell = \int_{t_0}^{t_1}\frac{v(t)\,\text{d}t}{a(t)}, $$ while the corresponding proper distance is $D = a(t_1)D_\text{c}$. For more details regarding co-moving and proper distance, see this post.

All that's left needed is an expression for $v(t)$. This is simply the SR formula for the relativistic rocket with constant proper acceleration $g\,$: $$ v(t) = \frac{g(t-t_0) + w_0}{\sqrt{1+[g(t-t_0) + w_0]^2/c^2}}, $$ where $$ w_0 = \frac{v_0}{\sqrt{1-v_0^2/c^2}}, $$ and $v_0$ is the initial velocity at time $t_0$; see this post for the derivation. By inserting the formulae for $a(t)$ and $v(t)$ in the integral above, we can calculate the travelled co-moving distance $D_\text{c}$. Also, the proper time elapsed on board is $$ \tau = \int_{\tau_0}^{\tau_1}\text{d}\tau = \int_{t_0}^{t_1}\sqrt{1-v(t)^2/c^2}\,\text{d}t. $$

Someone recently asked a question about numerically integrating @Pulsar's equations.

I don't agree with @Pulsar's answer from 6 years ago. You can't use the flat-space formula for $v(t)$; you need to derive how the rocket moves under constant proper acceleration in a Friedmann universe.

Take the spatially-flat Friedmann metric to be

$$d\tau^2=dt^2-a(t)^2(dx^2+dy^2+dz^2)\tag{1}$$

(with $c=1$). We will need the following Christoffel symbols:

$$\Gamma^t_{tt}=0$$

$$\Gamma^t_{tx}=0$$

$$\Gamma^t_{xx}=a\dot{a}\tag{2}$$

$$\Gamma^x_{tt}=0$$

$$\Gamma^x_{tx}=\dot{a}/a$$

$$\Gamma^x_{xx}=0$$

An overdot will consistently mean $d/dt$, not $d/d\tau$. We have to parameterize everything with $t$ because of the scale factor $a(t)$.

Consider a rocket that takes off at $t=t_0$ and accelerates from rest with constant proper acceleration $g$ in the $x$-direction. What is its coordinate velocity $v(t)=dx(t)/dt$?

Along its worldline, we have

$$d\tau^2=dt^2-a(t)^2dx^2\tag{3}$$

or

$$d\tau=dt\sqrt{1-a(t)^2v(t)^2}\tag{4}$$.

Introduce

$$\gamma(t)\equiv\frac{1}{\sqrt{1-a(t)^2v(t)^2}}\tag{5}$$

so that

$$v(t)=\frac{\sqrt{\gamma(t)^2-1}}{\gamma(t)\,a(t)}.\tag{6}$$

we will find that there is a relatively simple equation for $\gamma(t)$.

The two nonzero components of the rocket's four-velocity are

$$u^t\equiv\frac{dt}{d\tau}=\gamma\tag{7a}$$

$$u^x\equiv\frac{dx}{d\tau}=\frac{dt}{d\tau}\frac{dx}{dt}=\gamma v=\frac{\sqrt{\gamma^2-1}}{a}\tag{7b}$$

We double-check that

$$u\cdot u=(u^t)^2-a^2(u^x)^2=1.\tag{8}$$

The four-acceleration is given by

$$A^\lambda=\frac{Du^\lambda}{d\tau}=\frac{du^\lambda}{d\tau}+\Gamma^\lambda_{\mu\nu}u^\mu u^\nu.\tag{9}$$

Its two nonzero components are

$$\begin{align} A^t&=\frac{du^t}{d\tau}+\Gamma^t_{tt}(u^t)^2+2\Gamma^t_{tx}u^tu^x+\Gamma^t_{xx}(u^x)^2\\ &=\frac{dt}{d\tau}\frac{du^t}{dt}+a\dot{a}(u^x)^2\\ &=\gamma\dot{\gamma}+H(\gamma^2-1) \end{align}\tag{10a}$$

$$\begin{align} A^x&=\frac{du^x}{d\tau}+\Gamma^x_{tt}(u^t)^2+2\Gamma^x_{tx}u^tu^x+\Gamma^x_{xx}(u^x)^2\\ &=\frac{dt}{d\tau}\frac{du^x}{dt}+\frac{\dot{a}}{a}u^tu^x\\ &=\gamma\frac{d}{dt}\left(\frac{\sqrt{\gamma^2-1}}{a}\right)+\frac{\dot{a}}{a}\gamma\frac{\sqrt{\gamma^2-1}}{a}\\ &=\frac{1}{a}\left(\frac{\gamma^2\dot{\gamma}}{\sqrt{\gamma^2-1}}+H\gamma\sqrt{\gamma^2-1}\right) \end{align}\tag{10b}$$

where

$$H(t)\equiv\frac{\dot{a}(t)}{a(t)}.\tag{11}$$

Proper acceleration is the magnitude of the four-acceleration, so if the rocket has proper acceleration $g$, then

$$A\cdot A=(A^t)^2-a^2(A^x)^2=-g^2.\tag{12}$$

After some algebra, this reduces to the equation

$$\frac{\gamma\dot{\gamma}}{\sqrt{\gamma^2-1}}+H\sqrt{\gamma^2-1}=g\tag{13}$$

or, restoring $c$,

$$\gamma\dot{\gamma}=\frac{g}{c}\sqrt{\gamma^2-1}-H(\gamma^2-1).\tag{14}$$

In conjunction with the initial condition $\gamma(t_0)=1$, this first-order differential equation determines $\gamma(t)$. Note that $H$ is a function of $t$ which depends on the cosmological model one uses for the universe.

Once one specifies $a(t)$, and thus $H(t)$, one can solve this equation (presumably numerically) for $\gamma(t)$, and then find the coordinate velocity of the rocket,

$$v(t)=\frac{\sqrt{\gamma(t)^2-1}}{\gamma(t)\,a(t)},\tag{15}$$

the proper time $\tau(t)$ that elapses on the rocket,

$$\tau(t)=\int_{t_0}^t\frac{dt}{\gamma(t)},\tag{16}$$

and the comoving distance and proper distance that Pulsar mentions.

In flat space, $H=0$, so the equation I derived reduces to

$$\gamma\dot{\gamma}=\frac{g}{c}\sqrt{\gamma^2-1}.\tag{17}$$

Taking $t_0=0$, it has the solution that Baez mentions,

$$\gamma=\sqrt{1+(gt/c)^2}.\tag{18}$$

Unfortunately, it also has the trivial solution $\gamma=1$. There are two solutions to this first-order ODE with the same initial condition! This can happen with first-order differential equations of the form $dy/dx=F(x,y)$ when $F$ or $\partial F/\partial y$ is not continuous, as it is in this equation. (This lack of a unique solution is discussed here.)

This causes problems when trying to solve the equation with $H$ numerically: It also has the trivial solution $\gamma=1$, and Mathematica's NDSolve[] just finds this trivial solution.

The fix is to change variables from $\gamma$ to $\lambda$ where

$$\lambda=\sqrt{\gamma^2-1}.\tag{19}$$

The differential equation then takes the remarkably simple and well-behaved form

$$\dot{\lambda}(t)=\frac{g}{c}-H(t)\lambda(t).\tag{20}$$

The flat-space solution is obviously $\lambda(t)=gt/c$. It is easy to use Mathematica's NDSolve[] to solve it in an expanding universe with some $H(t)$. After solving for $\lambda(t)$ with the initial condition $\lambda(t_0)=0$, use

$$\gamma(t)=\sqrt{1+\lambda(t)^2}\tag{21}$$

to get $v(t)$, $\tau(t)$, etc. as discussed above.

The cosmological model typically used today is the Lambda-CDM model. In this model, the current density of dark energy as a fraction of the critical density is $\Omega_\Lambda=0.6911$ and the current fraction for matter (both dark and non-dark) is $\Omega_m=0.3089$. Since these add up to 1, the fractions for radiation and curvature can be taken to be zero: $\Omega_r=\Omega_c=0$. There turns out to be an analytic solution of the Friedmann equation:

$$a(t)=\left(\frac{1-\Omega_\Lambda}{\Omega_\Lambda}\right)^{1/3}\sinh^{2/3}\frac{t}{t_\Lambda}\tag{22}$$

where

$$t_\Lambda=\frac{2}{3H_0\sqrt{\Omega_\Lambda}}.\tag{23}$$

This gives

$$H(t)=\frac{\dot{a}(t)}{a(t)}=\frac{2}{3t_\Lambda}\coth{\frac{t}{t_\Lambda}}.\tag{24}$$

Thus, for numerical calculations, the only other parameter besides $\Omega_\Lambda=0.6911$ that one needs is the Hubble constant today. There is some dispute about its value, but the Wikipedia article on Lambda-CDM takes it to be

$$H_0=67.74\,\text{km}\,\text{s}^{-1}\,\text{Mpc}^{-1}.\tag{25}$$

This corresponds to

$$1/H_0=14.44\,\text{Gy}\tag{26}$$

and gives

$$t_\Lambda=11.58\,\text{Gy}.\tag{27}$$

I was unable to find a reference to a similar analysis, so I worked this out myself. The most interesting consequence of my analysis is that

$$\lim_{t\to\infty}\lambda(t)=\lim_{t\to\infty}\frac{g}{cH(t)}=\frac{g}{cH_0\sqrt{\Omega_\Lambda}}=\frac{3g\,t_\Lambda}{2c}\approx1.79\times 10^{10}.\tag{28}$$

(Just set $\dot{\lambda}=0$.)

So, in the Lambda-CDM universe, an accelerating-at-$1g$ rocket's $\gamma$ does not approach $\infty$ but rather a very large number.

If any reader is aware of a similar analysis, I would appreciate a reference to it.