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?