What does an outgoing null geodesic for FLRW metric look like? I know it involves setting the line element to zero, as well as setting $d\phi$ and $d\theta$ to zero. But other than that I do not know what to do.
I get an expression involving $a(t)$, the scale factor. I do not know how to integrate it.
 A: The metric is, if we set the angles to zero, $$ds^2 = c^2 dt^2 - a(t)^2 d\chi^2$$ where I use co-moving coordinates $\chi$. Setting $ds^2=0$ just gives the differential equation $$d\chi^2 = \frac{c^2}{a(t)^2}dt^2.$$ Taking the square root, $$d\chi = \pm   \frac{c}{a(t)}dt.$$ Now we can integrate and get $$\int d\chi=\int \pm   \frac{c}{a(t)}dt$$ which implies $$\chi(t) = \pm c\int_0^t    \frac{dt}{a(t)}$$ (there was a constant of integration, but since it makes sense to set $\chi(0)=0$ it vanishes).
This is as far as we can get without more information about $a(t)$. Unfortunately for realistic models it does not have a closed form. The good news is that we have computers that allow us to solve this numerically. The bad news is that the equation for $a'(t)$ is somewhat unruly. The good news to that is that one can use an implicit formula to get values of the time $t$ when $a(t)$ had a given value and then use them for integrating $\chi(t)$ numerically:
$$t(a) = H_0 \int_0^a \frac{d\alpha}{\sqrt{\Omega_m \alpha^{-1} + \Omega_k + \Omega_\Lambda \alpha^2}}$$
The bad news with that formula is that it is still not a closed 
formula (unless you like elliptic functions or use a very special case).
A nice special case is when only $\Omega_\Lambda$ matters, giving $t(a)\propto \log(a)$ or $a(t)=e^{kt}$ (a flat dark energy dominated universe; $k$ depends on the expansion rate). In this case the geodesic is simply $$\chi(t) = \pm c \int_0^t e^{-kt} dt =\pm \left(\frac{c}{k}\right) (1-e^{-kt}).$$ This shows that in co-moving coordinates outgoing geodesics will never pass beyond a certain distance. In proper distance coordinates $x(t)=a(t)\chi(t)$ they of course head off towards infinite distance due to the expansion.
