Momentum of a moving object in FRW metric according to an observer comoving with cosmic expansion I would like to show that in an FRW metric the momentum of a freely falling object decays as the inverse of the scale factor. I know there are many proofs and arguments for this but I am trying to get this starting from geodesics and having some trouble.
General logic
My general approach is the following.
The object is moving along some path. I want to find the momentum of the object that would be measured by a coincident, locally inertial observer, i.e. a detector that is freely falling and locally cannot tell that spacetime is curved.
In the locally inertial frame where the observer is momentarily at rest the 4-velocity of the observer has components $u^\mu\equiv d\xi^\mu/d\sigma= (1,0,0,0)$, where $\xi^\mu$ are the coordinates of this local inertial (cartesian) frame and $\sigma$ is the observer's proper time (in this frame $d\sigma=d\xi^0$). In this same frame the object has 4-momentum $p^\mu \equiv md\xi^\mu/d\tau$, where $\tau$ is the object's proper time, and $m$ the rest mass.
Since special relativity holds in this frame the components of $p^\mu$ have the values $p^\mu=(E_\mathrm{obs},\mathbf{p}_\mathrm{obs})$, the energy and momentum that would be measured by the observer at rest in this frame. They obey $E_\mathrm{obs}^2-\mathbf{p}_\mathrm{obs}^2=m^2$. We can retrieve the energy $E$ by contracting $u^\mu$ with $p^\mu$:
$$
E_\mathrm{obs}=g_{\mu\nu}u^\mu p^\nu,
$$
where $g_{\mu\nu}=\text{diag}(1,-1,-1,-1)$ are the components of the metric in this locally inertial coordinate system.
The magnitude of the 3-momentum in this system is just
$$
|\mathbf{p}_\mathrm{obs}| = \sqrt{(g_{\mu\nu}u^\mu p^\nu)^2 - m^2}.
$$
The right hand side is invariant under general coordinate transformations. Therefore, I think we can calculate the RHS in any coordinate system we want and the value will be the magnitude of the 3-momentum as measured by the observer  corresponding to the path $u^\mu$ in a locally inertial frame.
Is this ok so far?
FRW metric
$x^\mu = (t, \mathbf{x})$ are the coordinates of the FRW metric, defined by:
$$
d\tau^2 = dt^2 - a(t)^2 \left(d\mathbf{x}^2 + K \frac{(\mathbf{x}\cdot d\mathbf{x})^2}{1-K\mathbf{x}^2} \right),
$$ where $K =0$, $+1$, or $−1$.
The geodesic paths can be found using the Euler-Lagrange equations (see this post). The resulting equations are
\begin{align}
0 &= \frac{d^2 t}{d\tau^2} + a\dot{a} \left(\mathbf{x}'^2 +\frac{K(\mathbf{x}\cdot \mathbf{x}')^2}{1-K \mathbf{x}^2}\right)\\
\mathbf{0} &= \frac{d}{d\tau} \left[ a^2\left(\mathbf{x}' + \frac{K(\mathbf{x} \cdot \mathbf{x}')\mathbf{x}}{1-K\mathbf{x}^2}\right) \right] - \frac{K (\mathbf{x} \cdot \mathbf{x}')}{1-K\mathbf{x}^2} a^2 \left(\mathbf{x}' + \frac{K(\mathbf{x} \cdot \mathbf{x}')\mathbf{x}}{1-K\mathbf{x}^2}\right),
\end{align}
where a prime means $d/d\tau$ and $\dot{a}=da/dt$ is a function of $t$.
The simple solutions to these equations are $t=\tau, \mathbf{x}=\text{const}$. These solutions correspond to the "comoving" observers that move along with the cosmic expansion. The world lines of these observers all have 4-velocity $u^\mu=(1,0,0,0)$ in this coordinate system. When the object is at some location $\mathbf{x}$ at time $t$ I want to get the momentum as measured in a locally inertial frame in which the comoving observer that sits at position $\mathbf{x}$ is momentarily at rest.
To get the equations of motion for an arbitrary freely falling object you have to integrate the above equations. I was able to do this only for the second equation but I think that's enough.
Define $\mathbf{f}$ by
$$
\mathbf{f} \equiv a^2\left(\mathbf{x}' + \frac{K(\mathbf{x} \cdot \mathbf{x}')\mathbf{x}}{1-K\mathbf{x}^2}\right),
$$
to write the second equation as
$$
\mathbf{0} = \frac{d \mathbf{f}}{d\tau} - \frac{K (\mathbf{x} \cdot \mathbf{x}')}{1-K\mathbf{x}^2} \mathbf{f}.
$$
Multiply through by the integrating factor $\sqrt{1-K\mathbf{x}^2}$ to get
\begin{align}
\mathbf{0} &= \sqrt{1-K\mathbf{x}^2}\frac{d \mathbf{f}}{d\tau} - \frac{K (\mathbf{x} \cdot \mathbf{x}')}{\sqrt{1-K\mathbf{x}^2}} \mathbf{f} \\
 &= \sqrt{1-K\mathbf{x}^2}\frac{d \mathbf{f}}{d\tau} + \mathbf{f} \frac{d}{d\tau}\sqrt{1-K\mathbf{x}^2} \\
&= \frac{d}{d\tau}\left(\mathbf{f} \sqrt{1-K\mathbf{x}^2}\right).
\end{align}
The solution is 
\begin{align}
\frac{\mathbf{c}}{\sqrt{1-K\mathbf{x}^2}} = a^2\left(\mathbf{x}' + \frac{K(\mathbf{x} \cdot \mathbf{x}')\mathbf{x}}{1-K\mathbf{x}^2}\right),
\end{align}
for some constant 3-vector $\mathbf{c}$.
The term $(\mathbf{x}\cdot\mathbf{x}')$ can be isolated by dotting this equation with $\mathbf{x}$:
\begin{align}
\frac{\mathbf{c}\cdot\mathbf{x}}{\sqrt{1-K\mathbf{x}^2}} &= a^2\left((\mathbf{x} \cdot \mathbf{x}') + \frac{K(\mathbf{x} \cdot \mathbf{x}')\mathbf{x}^2}{1-K\mathbf{x}^2}\right)\\
&= a^2 (\mathbf{x} \cdot \mathbf{x}') \left(1+ \frac{K\mathbf{x}^2}{1-K\mathbf{x}^2}\right)\\
&= a^2 \frac{\mathbf{x} \cdot \mathbf{x}'}{1-K\mathbf{x}^2}.
\end{align}
Now $\mathbf{x}'$ can be written in terms of $\mathbf{x}$:
\begin{align}
a^2 \mathbf{x}' &= \frac{\mathbf{c}}{\sqrt{1-K\mathbf{x}^2}} - a^2 \frac{K(\mathbf{x} \cdot \mathbf{x}')\mathbf{x}}{1-K\mathbf{x}^2}\\
&= \frac{\mathbf{c}}{\sqrt{1-K\mathbf{x}^2}} -
\frac{K(\mathbf{c}\cdot\mathbf{x})\mathbf{x}}{\sqrt{1-K\mathbf{x}^2}}\\
&= \frac{\mathbf{c} - K(\mathbf{c}\cdot\mathbf{x})\mathbf{x}}{\sqrt{1-K\mathbf{x}^2}}.
\end{align}
Finally the spatial part of the object's 4-momentum is
\begin{align}
p^i = m\frac{d\mathbf{x}}{d\tau} = m \frac{\mathbf{c} - K(\mathbf{c}\cdot\mathbf{x})\mathbf{x}}{a^2\sqrt{1-K\mathbf{x}^2}}.
\end{align}
Now to try to use the $|\mathbf{p}_\mathrm{obs}| = \sqrt{(g_{\mu\nu}u^\mu p^\nu)^2 - m^2}$ formula I had earlier.
Since $u^\mu=(1,0,0,0)$, $g_{00}=1$, and $g_{0i}=0$ at $\mathbf{x}$, $|\mathbf{p}_\mathrm{obs}| = \sqrt{(p^0)^2 - m^2}$.
Writing out $g_{\mu\nu}p^\mu p^\nu = m^2$ gives
\begin{align}
m^2 &= (p^0)^2 - m^2 a^2\left(\mathbf{x}'^2 + K \frac{(\mathbf{x}\cdot \mathbf{x}')^2}{1-K\mathbf{x}^2} \right) \\
&= (p^0)^2 - m^2 (\mathbf{f} \cdot \mathbf{x}') \\
&= (p^0)^2 - m^2 \frac{\mathbf{c}}{\sqrt{1-K\mathbf{x}^2}} \cdot \frac{\mathbf{c} - K(\mathbf{c}\cdot\mathbf{x})\mathbf{x}}{a^2\sqrt{1-K\mathbf{x}^2}} \\
& = (p^0)^2 - m^2 \frac{\mathbf{c}^2 - K(\mathbf{c}\cdot\mathbf{x})^2}{a^2 (1-K\mathbf{x}^2)},
\end{align}
and I get
$$
|\mathbf{p}_\mathrm{obs}| = \sqrt{(p^0)^2 - m^2} = \frac{m}{a}\sqrt{\frac{\mathbf{c}^2 - K(\mathbf{c}\cdot\mathbf{x})^2}{1-K\mathbf{x}^2}}.
$$
But I don't think this is proportional to $1/a$ in general.
If the object passes through the origin then everything works as I expect. In this case $\mathbf{c}$ is parallel to $\mathbf{x}'$ when $\mathbf{x}=0$ and the equation of motion keeps it moving along a straight line parallel to this initial $\mathbf{x}'$ so $\mathbf{x}$ remains parallel to $\mathbf{c}$ and $(\mathbf{c} \cdot \mathbf{x})^2 = \mathbf{c}^2\mathbf{x}^2$ and the term in the square root is just a constant $|\mathbf{c}|$.
But for an object flying somewhere else, not passing through the origin I think the angle between $\mathbf{c}$ and $\mathbf{x}$ will change as it goes (and in particular it will not always be $0$) and so the momentum observed by a locally inertial observer will not simply scale as $1/a$.
I guess maybe the angle between $\mathbf{c}$ and $\mathbf{x}$ changes as $\mathbf{x}^2$ changes in just such a way as to keep the term in the square root constant but I don't see how to show that.
Or do I have the entire logic of the exercise wrong?
 A: It turns out that last factor in the square root is actually a constant -- so everything works nicely and momentum decays as $1/a$ for any object moving on a geodesic as it's supposed to.
I took the derivative of the term in the square root:
$$
\frac{\mathbf{c}^2-K(\mathbf{c}\cdot \mathbf{x})^2}{1-K\mathbf{x}^2} \equiv A.
$$
This introduces the terms $\mathbf{c}\cdot \mathbf{x}'$ and $\mathbf{x}\cdot \mathbf{x}'$, which we can get rid of using some of the formulas in the original post:
\begin{align}
\mathbf{c}\cdot \mathbf{x}' &= \frac{\mathbf{c}^2-K(\mathbf{c}\cdot \mathbf{x})^2}{a^2\sqrt{1-K\mathbf{x}^2}} = \frac{\sqrt{1-K\mathbf{x}^2}}{a^2}A \\
\mathbf{x}\cdot \mathbf{x}' &= \frac{1}{a^2}\sqrt{1-K\mathbf{x}^2}(\mathbf{c}\cdot \mathbf{x})
\end{align}
\begin{align}
\frac{dA}{d\tau} &=\frac{-2K (\mathbf{c}\cdot \mathbf{x})(\mathbf{c}\cdot \mathbf{x}')}{1-K\mathbf{x}^2} + \frac{\left(\mathbf{c}^2-K(\mathbf{c}\cdot \mathbf{x})^2\right)2K(\mathbf{x}\cdot \mathbf{x}')}{(1-K\mathbf{x}^2)^2}\\
&= \frac{-2K (\mathbf{c}\cdot \mathbf{x})}{1-K\mathbf{x}^2}\frac{\sqrt{1-K\mathbf{x}^2}}{a^2}A + \frac{A}{1-K\mathbf{x}^2} \frac{2K}{a^2}\sqrt{1-K\mathbf{x}^2}(\mathbf{c}\cdot \mathbf{x})\\
&= 0.
\end{align}
So 
$$
|\mathbf{p}_\mathrm{obs}| = \frac{m}{a}\sqrt{A}= \frac{\text{const}}{a}.
$$
