Wave packet in curved spacetime It is known that the classical equation of motion for a scalar field wave packet on a curved spacetime background gives the geodesic trajectory (the e.o.m. is $(\nabla_\mu \nabla^\mu + m^2) \Phi=0$). However, I couldn't see that. 
How can one derived the geodesic equation from the above e.o.m. ?
 A: This is covered (in two pieces) in Chapter 4 of Wald's General Relativity for the Maxwell field;  I'll adapt his proof here.
We start with the equation $(\nabla^a \nabla_a + m^2) \Phi = 0$.  Let us look for a solution to this equation of the form
$$
\Phi = \Phi_0 e^{iS},
$$
where $S$ is a real-valued function and $\Phi_0$ is "slowly varying" compared to $S$.  This is the standard sort of "geometrical optics" approximation;  in practical terms,  it means that we can ignore any derivatives of $\Phi_0$.  Thus, we have
$$
\nabla_a \Phi = \Phi_0  e^{iS} (i \nabla_a S)\quad \Rightarrow \quad \nabla^a \nabla_a \Phi + m^2 \Phi = \Phi_0  e^{iS} \left[ m^2 - (\nabla^a S) (\nabla_a S) + i \nabla^a \nabla_a S \right] = 0.
$$
Both the imaginary and real parts of the quantity in square brackets above must vanish:
$$
m^2 = (\nabla^a S) (\nabla_a S), \qquad \nabla^a \nabla_a S = 0.
$$
If we let $k_a = \nabla_a S$, the first equation implies that $k_a k^a = m^2$, as we expect.  Moreover, if we take the derivative of both sides of the first equation, we get
\begin{align*}
0 &= \nabla_b \left[ (\nabla^a S) (\nabla_a S) \right] \\
&= 2 \nabla^a S \nabla_b \nabla_a S && \text{(product rule)} \\
&= 2 \nabla^a S \nabla_a \nabla_b S && \text{(derivatives commute on scalars)} \\
&= 2 k^a \nabla_a k_b,
\end{align*}
and so the vector $k_a$ satisfies the geodesic equation.
A: As in your comment, make an Ansatz $$ \phi = \exp(i f)$$ where $f$ is an unknown scalar function such that $(df)_\mu (df)^\mu = m^2$. Then the equation of motion for $\phi$ is equivalent to the geodesic equation for $(df)_\mu$. 
