 Solution

Divide the time component of the geodesic equation by $(P^0{)}^2$ and use $d/d\lambda =P^{0}d/d\eta$. To first order, $P^i/P^0={\hat{p}}^i(1+\delta N-\hat{\mathbf{p}}\cdot \nabla \psi )$. Substitution of the supplied Christoffel symbols makes the background $2\mathcal{H}$ terms cancel the derivative of $a^{-2}$. The derivatives of the lapse and shift also combine, leavingThus the comoving photon energy is conserved in the unperturbed spacetime and changes through the gradient of the scalar gravitational source.