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. ?
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.
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$.
