Constant acceleration over cosmological distances At constant acceleration in special relativity, the time differs for a stationary observer and the astronaut. see the following article for an in-depth explanation:
Relativistic Rocket
However, when large distances are involved, due to the expansion of the universe, the article says that general relativity equations will have to be used instead. So what are the general relativity equations that should be applied to a relativistic constant acceleration involving large distances?
 A: 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.
$$
A: 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.
