On the other side of the question, do you have any Riemannian manifold for which your conjecture does not hold?

 Dear  Prof.  Leyvraz   Thank  you  for  your  answer    After  your  answer, I  search  more and I  understand  from  proposition  1.8  of  this  file that the  solutions  are  always  geodesics(At  least  locally).

On  the other  hand  it  is  difficult  to  imagine  this  fact  for the  torus,  since  the  gradient  vector  field  has  singularities  of  negative  index, say  saddle  point. It  is  difficult  to imagine  that  the  geodesics  around  a  saddle  point  are so  divergent  from  each other

http://cedricvillani.org/wp-content/uploads/2012/08/P12.CIME_.pdf

The distance d^2(x,p) is ``multivalued'' in the large on manifolds like a torus, so any arguments in the large are tricky.

This property is true for any Riemannian manifold. If you use the function d instead of d^2, the trajectories of the gradient (whereever d is smooth) are geodesics parametrized by arc length.

Best regards

Jost

Hi, If \nabla denotes the Levi-Civita connection on he Riemannian manifold (M,g) and D/dt denotes the related covariantt derivative, and the trajectories x(t) of the function d^2(x,p) are geodesics (after a reparametrization) then in any local chart on M, (D/dt v(t)) is proportional to v(t), where, v(t)=d/dt(d(x(t)).  Moreover, the function f(x):=d^2(x,p) is a solution to the following PDE

                            \nabla_i \nabla_j f=\lambda(f) g_ij,

where, \lambda is an appropriate function defined on M and also is a function of f. This manifolds are studied upon special forms of \lambda (e.g polynomial on f). See for example

Y. Tashiro, Complete Riemannian manifolds and some vector fields, Trans. Amer. Math. Soc., 117 (1965), pp. 251–275.

The above study may help you to discover the Riemannian manifold satisfying your property (f may be replaced by d^2(x,p)).

All the best